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SYMBOLS USED 


A : 

Plate thickness 

(Nondimensional length unit) 

b : 

A half of the layer thickness 

N: 

Layer number in 

the plate 

p: 

Density 


t : 
o 

Impact time 



P (x,,t); Impact stress 
o i 

* x^(n)» X 2 (?): Coordinate variables 

* t(t): Time variable 


T : Nondimensional time unit = a/»x77/p 

o 6o 


* t(T), <^2.1* Stress tensor and its components 

* u^, u(U), v(V): Displacement vector and its components 

♦ Elastic Moduli 


'•ijU’. ij^'ij 
X, p: Lame's constants 


e . . , e.: Infinitesimal strain tensor 

ij 1 

; Laplace transform (in x) and Fourier transform (in n) of A 

P (5): Legendre polynomial of C 

n 

* k(<) : Wave number 

* aj(iu) : Frequency 

0, a, B: Phase shifts (wave number through the thickness) 

Cjj, Cgi Dilatational and shear wave speed 

G*(oj) = G'((jj) + iG"(a>); Complex modulus of elastomer 

D: Thickness of viscoelastic layer 

h: A half of the crack length 


Quantities in ( ) are nondimensional quantities. 


ii 




Preface and Summary 


This report is the last of a series on the response of composite 
plates to impact forces. The motivation for these studies has been an 
attempt to understand the damage to aircraft jet engine fan blades by 
foreign object impact such as ice balls, stones, and birds. In addition, 
the National Aeronautics and Space Administration, sponsors of this 
research, have sought to develop computer codes from these analyses which 
will aid the fan blade designer in locating potential failure modes and posi- 
tions and thus assist in optimizing fan blade fabrication to create greater 
impact tolerance. 

The basic approach of the principal investigator in these studies has 
been to use wave propagation techniques to model the early response of 
composite plates to impact type forces. In using the wave method, the 
plate can be simplified in the analyses by neglecting reflections from edge 

t 

boundaries far from the impact point. fhUB?,whlle”the’'ovefall-'-geometry ) 

^thematical->^iK^d^sTieari:^e "point '^jf^^^iinpa'ct^'liave ^b 

The basic model for the composite plate studies has been the anisotropic plate 
theory as extended by Mindlin [ 1 ] to account for wave phenomena. The 
plate equations were used as an approximation of the exact theory of 
elasticity because they lead to simpler computational schemes for evalua- 
ting average stresses and displacements in the plate. 

Fourier and Laplace transform techniques have been used throughout 
these studies and inversion of the transforms has been accomplished 
with a fast Fourier transform algorithm. This algorithm is an effective compu- 
tational tool but requires careful scaling of the impact problem in both space 
and time 
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variables. When it is properly used it can lead to calciilations of thousands of 
stress values in a fraction of the time of conventional finite element schemes. 

In summary, the use of plate models for the fan blade impact has avoided 
the analytical complexities of the exact theory of elasticity as well as 
the computational difficulties of finite element methods. 

In earlier reports both central and edge impact of an anisotropic 
plate were studied, [ 2 - U ]• In those reports only wave propagation in 
the plane of the plate was investigated. In another report [5] a 
multilayer plate model was developed in order to study impact induced 
wave propagation in both the thickness and in plane directions. In 
this final report further results are presented from the multilayer 
model. The composite plate has been modeled with as many as eight 
separate layers. Each layer may itself have several plys, so that 
effective anisotropic constants must be used for each layer in the analysis. 

The mathematical model exhibits wave propagation in both the thickness and 
inplane directions. Impact generated waves are shown to lead to inter- 
laminar shear stresses and Interlaminar tensile stresses during and after 
impact . 

This report also presents an analysis of an impact damping mechanism. 

This consists of thin damping layer introduced between composite layers in the 

mathematical model. The resulting response due to impact shows that considerable 
reduction of stress can be achieved. However it appears that this stress 
reduction is linked to the lower elastic moduli of the damping sublayers 
and not the viscous nature of the sublayer. 

Fianlly an attempt was made to analyze the impact of a plate with a 
crack. While the problem has been formulated, no progress was made on 
obtaining numerical answers to the crack problem. 
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I. INTRODUCTION 


The present research is a continuation of our previous work on the 
stress wave propagation in a laminated composite [2-5]. It has heen 
motivated by the problem of the impact on jet engine fan blades caused by 
ingestions of foreign materials, such as birds and hailstones. The 
successful application of fiber-reinforced composite materials depends on 
the ability of these materials to withstand forces due to such impact. 

The simplest approach to examine the dynamic behavior of a composite 
plate is based upon the work of White and Angona [ 6] . In their work, 
referred to as the effective modulus theory , the response of the composite 
plate to waves propagating normal to the laminate is predicted by a 
single constant wave speed, regardless of the internal structure of 
composites. Even though this theory is satisfactory for many problems, it 
fails in the case of some problems when the wave lengths become short. To 
overcome this limitation. Sun and et al. proposed a model which includes 
the effects of internal structure, such as the layer thickness [t 1- In their work, 
referred to so the effective stiffness theory , displacements of both the reinforc^ 
ing layer and the matrix layer are expressed as linear expansions about the 

midplanes of the layers and approximate equations of motion are derived 
for both layers. Then these approximate equations are required to satisfy 
the continuity conditions of displacement and stress components on every 
interface. Using this model the propagation of harmonic waves has been 
examined . 

More recently a number of researchers have presented models for multilayer 
plates either by the discrete-continuum theory or the continuum mixttire 
theory [8-lit]. Many, however examined only the frequency- 
wave number dispersion relationship and stopped short of the transient 
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impact problem except for a few experimental or numerical works using 
the finite element method which sometimes show a considerable discrepancy 
from the experimental results • 

In this report we present a new attempt to mathematically model the 
multilayer plate and develop a method by which we can examine the 
transient propagation of an impact wave in the plate, not only along 
the longitudinal direction but also through the thickness direction of the 
plate as well, using an inexpensive Fast Fourier Transform algorithm 

[3,15]. 

The composite plate under consideration for the first part of the 
present report is imagined to comprise N identical elastic layers. And 
each layer is made of a number of unidirectional plys lying alternately at 
a layup angle i<J> from the 'symmetry axis, as shown in Fig. 1. Then the 
elastic properties of the plate depend on the layup angle 4>. A key 
assumption for the first step of the work is that all the layers are 
identical. While restricting the application, this assumption allows us 
to formulate the problem using difference-differential equations due to 
a rather simple periodic structure of the plate. This technique for 
periodic structures has been widely used in the study of electrical 
transmission lines [I6] and in the vibration of multistory buildings [IT]. 
By means of an approximate plate theory of Mindlin [I6], a set of approx- 
imate equations of motion is developed for a typical layer using the inter- 
laminar stresses as explicit variables. The relative motion of a layer 
to the adjacent layers is related by phase shifts which represent the 
solution of the difference parts of equations. In this way the number of 
the layers can be increased without increasing the size of matrix in the 
numerical process of invert to satisfy the boundary conditions. 
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It is also well understood that a thin viscoelastic layer present 
between elastic layers can reduce the elastic impact energy significantly 
by dissipating the strain energy into heat [19,20], In our previous work 
[ 5 ] an elastomer layer is presented between a composite half space and a 
protection strip on the edge on which the impact is applied. Numerical 
results of the work showed an appreciable reduction in the normal stress. 

As an extension of this research and the first part of this report 
we now examine the wave propagation in a composite plate made of two 
elastic layers and an elastomer layer. Generalization of this problem is 
straightforward by assuming that our new periodic composite layer is now 
made of an elastic sublayer and a viscoelastic sublayer lying alternately. 

We can now develop the approximate theory which includes both sub- 
layers. For the second part of the present research we will examine the 
simplest case of this kind, i.e., an impact on a composite plate 
consisting of two elastic layers and an elastomer layer between them. 

Another possible extension of the multilayer composite plate which 
can be found in frequent practice is discussed in the last part of this 
report. In this chapter a crack is introduced on the interface 
between two elastic layers which represent the final step before a 
failure occurs in the composite either by spalling or by excessive shear 
stress. Such crack problems constitute the main part of the study of 
fracture mechanics, A serious mathematical difficulty arises even in the 
static problems because of the mixed boundary conditions along the crack 
direction. The difficulty becomes more serious in the case of dynamic problems 
due to the diffraction of waves at the crack tip [21-21+]. By 
employing the approximate equations of motion developed in the first 
part, the transient wave problem has been formulated and dual integral 
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equations are obtained after application of the mixed boundary conditions. 

But the resulting dual integral equations are not easy to solve and are under 
investigation at this time. 

In the results presented in this report only a line impact has been 
examined. This has simplified the calculations and saved computer time in 
testing the model. The technique, however, can be extended to the two- 
dimensional or central impact problem. Since the impact speed is very high 
(-450 m/sec) and the Impact time is short (< 100 ysec), the impact can be 
in the range of the elastic-plastic impact or even in the range of the 
hydraulic impact. But the initial transmission of impact energy is propa- 
gated by elastic waves, as if in anunbounded plate, and it is useful to 
investigate the problem by means of the linear theory of anisotropic 
elasticity in an infinite composite plate. 



II. IMPACT ON MULTILAYER ELASTIC PLATE 


1. Formulation 

Basic Theory of Linear Anisotropic Elasticity 

Cauchy's equations of motion in cartesian tensor form without body 
forces are given by 


.i 


pu. 


o . . = cr. , 

13 31 


(II-l) 


where the repeated index implies summation on that index. A comma repre- 
sents a partial differentiation with respect to the index following the 
comma and a superposed dot denotes a time derivative. 

tensor is related to the infinitesimal strain tensor e. . by 

i3 




or 


c . . E . 

i3 3 


(II-2) 


The condensed elastic moduli 
materials 



has the following form for orthotropic 


c . . 

13 


» ^12 * *”13 ’ ® ® ® 

* ^22 ' ^23 * ^ ^ ^ 

^13 ’ *^23 ’ *^33 ’ ° ° ° 

0 , 0 , 0 , *"44 > ® ® 

0 , 0 , 0 , 0 , Cj^ , 0 

0 , 0 , 0 , 0 , 0 , 
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Analysis of a Layer 

For a layer shown in Fig. 1 we employ the approximate plate theory 
of Mindlin (18 ] and the displacement field u is expanded in terms of 
the Legendre polynomials as 


u(x^,X2,X2,t) 


E u^“^(x ,x ,t) 
n=0 ~ 




(II-3) 


where 5 is the local coordinate along "the thickness direction, normalized 
by b, a half layer thick. 

Instead of solving Eq. (II-l) directly we solve a new approximate 
equation of motion which is obtained through a variational process by 
integration of Eq. (II-l) over the thickness 5 (see fl], [23]). The result is 


b 


aja 


+ [P„(C).a2^ 


] 


1 

C*-l 



2pb ..(n) . j = 1,2,3 
2n+l j ’ a = 1,3 


(II-A) 


where 


= f P (C)-o . dC 
aj n'^-' aj 


*(n) 


/ 

-1 


dP(C) 


2j dc ''2j 


dC 


By substituting the constitutive relation (II-2) for the displacement 
expansion (II-3) into the above approximate equations of motion, we can 
find governing equations of motion in terms of * • • 

The accuracy of this approximate theory depends on how many terms of the 
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displacement field we retain. Since the complexity in formulation 

increases rapidly with the number of terms included we keep terms only up 

to second order. Furthermore, we will examine harmonic waves propagating 

along the x^ and X 2 directions so that we drop terms and set 

-r — { } = 0. Next to get rid of the undesired coupling with higher modes 

3x^ 

we set ii^ = U 2 = 0. Then the resulting equations are 




2bpilf) 


^*’‘^66^b“l,l‘*^2,ll^ ^®22"®22^ 


= 2bpu 


( 0 ) 


f‘^^l“l!il'^12'^2!i^ ■ ^°2lSl^ =|bpu{^^ 


( 1 ) 


2b. (1) ^3 (2). (0)^1 (1), ^ N 2, ..(1) 

3 ^^66'^2,11 b^66^1,l^ ‘ ^^^12'^l,l‘^b'^22'^2 ^ ^'^22'^22^ - 3bpU2 


^®22"*^22^ 




^<‘=12“ia+b'=22“f^> 


(1)_^3 (2). 

^^‘^66^2,1'^b *^66“l ^ 


= 0 


= 0 


(II-5) 


where the sign + and - represent the stress components on the top and 
bottom surfaces of the layer under examination, i.e., at 5 = +1. Here 
we notice that the first, fourth, and last equations are written in terms 
of 5 where (n+i) is an odd integer and represents the thickness 

stretching motion (or symmetric motion) . In the rest of the equations in 
which (n+i) is an even integer the displacements represent the flexual 
motion (or antisymmetric motion). Hence, this process has decoupled the 
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stretching motion from bending motion I To get rid of the 2nd order 


( 2 ) 


and 


modes from Eq. (II-5) we solve the last two equations for 

( 2 ) 

, and insert them into the remaining equations. Then Eq. CII~5) can 
be reduced as follows: 




+ <“22- 


^22^ = 2bpU2°^ 


(II-6) 


2^ ..(1) O. r 1 

T^ll l,r^ 66^ b 


3^<“t2-‘’22’>l <“nSl' 


2u -(1) 

3bpui 


( 1 ) 

2 


-2<“i2“”M‘=22“2"> <“22S2' ‘ 1“’" 

where *” ^ii” ^12^^22 • 

Plate Analysis 

In view of the Legendre polynomial expansion, the displacements on 

the both sides of a layer can be written as uT = uf^^ ± uf^^ since the 

i X 1 

governing equations for a layer, Eq. (II-6), only include terms up to 
the first order of expansion, i.e., a linear expansion. Remembering that 
this analysis is valid for any arbitrary layer in a plate, say the 
nth layer, equation (II-6) can be immediately written as 


4r 

These two motions are. of course»coupled through the boundary conditions. 
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p (ii +ii . ) 

n n-i 


C, , (u +U , ) +~^ (v -V T (t -T ) 

11 n n-1 ,11 b n n-1 ,i b n n-1 


n n — 1 


-^c,-(u +u ,) --^c^_(v -V t)+‘t(o +o ,)+(t -t ) 
b 12' n n-1 ^ ^2 22 n n-1 b ' n n-1 n n-1'^ 


P(ii -ii -J 
n n-1 


c (u -u (u -u ,) --r c,,(v +v 

11 n n-1 ^2. ® n-1 b bo n n-1 ^ 


(II-7) 


+ — ^ (a -a ) + -^(t +t , ) 

C22 ^ jl ^ ° 


p (v_+v ) = c , , {-^(u -u - ) +(v +v ) } + -r(° 1 ) 

n n-1 66 b n n-1 n n— 1^11 d n n-1 


where o and t are used to represent 0^2 and 0^2 ^nd u and v 
denote and U 2 > respectively. These equations are the approximate 

equations of motion of a layer written in th e form of a difference-differential 
equation . For a plate made of N layer, the above equations contain 
4(N+1) unknowns (u^,v^,T^,a^, . . .Ujj,Vj^,Tj^,Oj^) and offer 4N equations. Since 
the additional four conditions are supplied by boundary conditions on 
the top and bottom surfaces, solutions of these equations can be found. 

In Eq. (II-7) we notice some important points. The first point is 
that the logitudinal coordinate and time variable are> continuous 

variables while the thickness coordinate X 2 is now discrete. This 
enables us to use integral transforms in x^ and time variables so that 
we can arrive at pure difference equations after integral transforms. 

The second point concerns the continuity conditions of stress and dis- 
placement. We note that u, v, -^22 ’ continuous across 

the layer boundary and these conditions are identically satisfied by 
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Eq. (II-7)- But the normal stress tangential to the layer boundary is not nec 

essarily continuous and Eq. (II-7) allows such a possibility. One can 
retain higher order terms in the displacement expansion given by Eq. 

(II-3) to give more accurate results. This can be achieved more easily 
by using Eq. (II-7) and increasing the number of layers in a plate under 
consideration. This process does not give any additional difficulties 
except a little more computer time. 


2. Dispersion Relationships of Harmonic Waves 
Harmonic Waves 

Before we examine the transient propagation of stress wave due to an 
impact we first investigate dispersion relations of harmonic waves in a 
composite plate governed by approximate equations of motion (II-7) . For 
harmonic waves propagating along the axis we assume 


{u ,v ,a ,T } 
n n n n 


{U ,V ,E ,T } e 
n n n n 


i(kXj^-o)t) 


(II-8) 


Substituting this into the approximate equations of motion (II-7) we obtain 


(m^-c K^)(U -KJ )+c iK(V -V )+b(T -T ,) = 0 

li n n~l 12 n n— 1 n n~l 


-3c ix(U +U )+(m^-3c„)(V -V„ )+ixb(T^-T^ ^)+3b(Z +E ,) = 0 

12 n n— 1 22 n n— 1 n n— 1 n n— 1 




(n-9) 




12 

+ iKb(I -E ,) = 

n n-1 


0 



11 


for n = 1, 2 . . .N. Here we set 


K = bk = k(-^), A = 2bN 


-2 ^2 2 2, A 2 

0) = pb 0) = p(jj (-^) 


and A is the total thickness of the plate. For a plate consisting of 

N layers, the boundary conditions require traction free surfaces, namely, 

T = Z = T„ = E„ = 0. When these conditions are applied to equation 

(II-9) we obtain 4N equations in terms of 4N unknowns 

U ,V ,T ,Z with n = 1,...N-1; U„,V„) . By setting the coefficient matrix 
n n n n N N 

to be singular, required dispersion relationships can be obtained. 


One-layer Plate 

The dispersion relationship for a plate made of a single layer can 

be found by setting N = 1 in equation (II-9) with E = T = E, = T, = 0. 

o o 1 1 

The resulting equations are now written in matrix form as follows: 


■ -2 2 
(o) ) 


C12IK 


-C12I-K ^(- 3 c.,.,+a!^) 


'22 


0 

0 


U, 

+ 

u 


0 




1 


0 



0 

0 


V. 


V 


0 




1 


0 





• 




— 



2 -2 







'^66^'' 

C66>c -0) 


^1 


u 

0 


0 

- 2 -2 








c,-ic +3c,,-(jj 

3c, ,iic 


V, 

+ 

V 


0 

11 66 

66 


_ 1 


0 




(II-IO) 


Then by setting the determinant of the coefficient matrix zero we obtain 
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21-2 -2 2 

CiiK^ - y(o) -3C22) (u) = 0 

2 1.-2 ^ 2 , .,-2 2 . - 

C66< - 3(c«> -c^i< - 3 ^ 65 ) -CggK ) = 0 


(II-ll) 


Here we notice that the first relationship corresponds to the state of 
deformation of U, = U and V, = -V , which represents the thickness 
extension of the plate (or the symmetric mode) , and the second describes 
the flexual deformation (or antisymmetric mode) . The exact theory of 
plates gives an infinite number of dispersion relationships, but because this 
model only has two inertia points (namely n = 0, l) , each of them having 
two components of displacement, we only have the first four relationships. 

Dispersion relationships and corresponding phase velocities for an 
isotropic plate with Poisson’s ratio 1/4 (namely X = y) are given in 
Fig. 2a and 2b up to the range where the wave length becomes equal to the 
plate thickness. Solid lines represent the symmetric modes and dotted 
lines the antisymmetric modes. As predicted by Mindlin and Medick the 
optical branch of the symmetric mode approaches the dilatation wave [I 8 ] . 

The acoustic branch of the antisymmetric mode starts from the bending 

motion and approaches the shear wave when the wave number k becomes 

* 

larger and larger . Similar relationships for an anisotropic plate made 
of 55% graphite fiber-epoxy matrix with a layup angle of 45° are shown in 
Fig. 3a and 3b. 


* 


See section 5 for discussion about the large wave nmber limit. 
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Two-layer Plate 

In this case we obtain eight equations by putting n = 1 and 2 in 
equation (II-9) . Boundary conditions require T = E = T_ = E„ = 0. By 

O Q ^ ^ 

following the same procedure we find the dispersion relations as 

r ,-2 2 , ,-2 „ , , 2 2,, -2 - 2 . .‘^ 66^^12 2 . 

{(m -c^^K )(o) -3 c22)-3c^2'" ^ 

-2 2 -2 2 -2 2 2 2 

+3(0) -C^^IC ){(0) -CggK )C0) “3Cgg)-3CgglC } 

,,-2 2,, -2 - 2 - . _ 2 2 , 

+3(o)^-CggK^) { (o)^-c^j^<^) (o)^-3c22)”3c^2'^^^ “ 0 

Again the first equation represents the symmetric mode and is shown as 
solid lines in Fig. 4 and 5. The second equation is plotted with dotted 
lines representing the antisymmetric mode. 

As expected we have six relationships since the this two-layer model is equiva- 
lent to a three-mass system with two degrees of freedom for each mass, 'i'/hen the 
wave number kA increases the following are observed: for the 

symmetric mode the upper optical branch approaches the dilatation wave^ whereas 

* 

for the antisymmetric mode the lower optical branch approaches the shear wave . 


= 0 
( 11 - 12 ) 


* See section 5 for discussions about the large wave nvunber limit. 



N-Layer Plate 


_2 

In general, we can obtain a 2(N+1) order poljmomial of o) by expand- 
ing a (UN) X (^+N) determinant and finding 2(N+1) dispersion relationships. 

Butj unfortunately, this process involves considerably complicated algebra 
and it may be necessary to develop a computer technique to find roots of 
an equation in determinant form (not in pol 3 nciomial form) . 

A difference equation approach can be used to solve the N set of 
four simultaneous first order difference equatiors given by Eq. (II-9) . 

This proceedure is neat and can be generalized for any number of layers 
as discussed in the next section; but the last step of this approach, where 
a long polynomial is to be solved again, is not any simpler than the previous 
direct method. 
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3. Impact on an Elastic Composite Plate 

Normalization and Integral Transforms of Governing Equations 

The governing equations given by (II-7) are first nondimensionalized 
as follows: 


{U ,V n) 
n n 


{u /A,U /A,x,/A} 
n n 1 


{C..,T ,E } = {c../c,,,T /c,,,o /c--} 

ij n n ij 66 n 66 n 66 


T = t/T 


where A is the total thickness of the plate and is the time required for 

the quasi-shear wave to travel the impact radius. Next we apply a Laplace 
transfom in x and a Fourier transform in n» i.e., 

00 

g(s) = / g(t)e ®^dx 

o 

00 

g(k) = f g(n)e^^'^dn 

/Ztt — « 


Then the resulting equations are 




C, 2 ik(U +U )-(^+2NC„)(V -V ,)+(!+! = 0 

iz n n-1 J IZ n-1 n-1 n n-1 6 N n n-1 


(n-13) 


-C ik(u -U T)-(fs +-:^k )(V+V -)+(E -E ,) = 0 

66 n n-1 ' 2N ' ' n n-1 n n-1 


-(H-+ +2NC,,)(U -U ,)+C ,ik(V+V ■) - - -(I -E ,)+(T -T ,) = 0 

3 6 N 66 n n-1 66 n n-1 6 NC 22 ^ 
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where the normalization factor f is given as 


f - J__A_ 

66 o 


bAp 

c 

66 o 


Solution of Difference Equations 

Since the simultaneous difference equations given are linear 
and all the coefficients are constants the solution [26] has to be 


{U ,V ,T ,E } 
n n n n 


= {A,B,C,D}e 


2i0n 


(11-14) 


where the phase shift 0 is complex, in general, and propagation through the 

thickness direction in the plate is characterized by 0. Namely, 0 is the 

wave number in the thickness direction. By substituting solution (11-14) 

into the difference equation (11-13) we obtain a set of four simultaneous 

homogeneous equations through which the relationships among the constants A,B,C, 

and D have to be determined. If we set the resulting coefficient matrix of 
A,B,C, and D to be singular we obtain the following equation for phase shift 0: 


cos^(fs^+^^)(fs2+^ k^ 


+ sin*^0('^® 


+ 2NC -^).(lii+^kW - k2) 

3 ^22 6 N 3 6 N ^^^^^66 6 NC 22 ^ ^ 




«V)} 


. 2 , 2 


2 . ^66, 2, ,fs^ . ^11, 2 




6N 


+2NC 


2 2 
cj k 

+ Ti^)] 


66 6 NC 


22 


(11-15) 


= a. 


4 

cos 0 


2„ 

a2Cos 0 


+ a. 


0 . 
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This equation implies that for a given wave number k along and a 

frequency s s represents the frequency for the case of harmonic waves), 
an infinite value of wave numbers exists for propagation through the 
thickness direction, but onlyfour of them are sufficient to give all linearly 
independent solutions of the form of Eq. (11-14). If we denote the solution 
of the phase shift equation as 


-a„+v^?-4 


-a2+/a2-4a^a3 

cos 6 = 


Jl. 

cos a = T 

2a, 


(11-16) 


the complete general solutions of difference equation (11-13) are 




— ~ 




— — 


- 

u 


A, 


Ao 




A, 

n 


1 


2 


3 


4 

V 


B, 


B, 


B~ 


B, 

n 


1 

2iBn 

i. 

-2iBn_^ 

3 

2 ion. 

4 




e + 


e + 


e + 


T 


C, 






c. 

n 


1 


2 


3 


4 

Z 


D, 


D., 


D., 


D, 

n 


1 


2 


3 


4 

- 


- 


— — 


_ 


^ __ 


-2ian 


(H-17) 

Next, by substituting the above solutions into the original difference 


equations (11-13) we find the relationships among A^, C^, and 


The results are 



Mi> Hi> <Ii> d> 
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X^( 3 ) 


X 2 ( 6 ) E 2 


X 3 ( 3 ) E 3 


• cos2nB+i 


X^CB) E 2 


X 2 ( 6 ) E^ 


X3(B) E3 


- sln2nB 


(11-18) 


Y^(a) E^ 


Y^Ca) E 3 


Y 2 (a) E 3 


Y 2 (cc) E^ 


• cos2na+i 


sln2na 


^3(0) E^ 


Y 3 (o) E 3 


where E^^ = + D2, E2 * - D2» E3 = C3 + C^, E^ *» C3 - and 


X^(B) = 


A^(B) 


(fs^+ -^^)(fs^ + -^k^)cos^B 
zn zn 


C 2 C 

+ { (if-+-^W 5j)-Cg5Cj2k2-c2jk2 jsln^B-cosB 


c c 


(11-19) 


A 2 ( 3 ) 


A 3 (B) 


. . 3„r^66^12, 2 . ^11, 2,„ „ ,, 2„ . 2.^11, 2, 

= 1 sxn k -(~3 — '■"5^1^ +2nCgg)}-cos BsxnB(fs + k ) 
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and 


A, (a) 

Y (a) 

A (a) 


A(a) 


2 2 2 
i cos a*sina{Cggk 





2 , 2 


2”<=66+6S 


k 


-)} + i sin^a(-^ 


22 


+ 2nC22) 


,<=66<^12'‘' fs! fulL , , , 

^ ^ 6 nC 22 ~ 3 6 n 


Aj^(a) = 


3 2 

cos a(fs +■ 


166 2 

2 n ^ ^ 


+ sin2a.cosa{-(^)^^(fs^+-^k2) +^c^^+(i±-+ 2nc^^)) 


A2<a) 


k-cos asina(C^ 2 '*^gg) 


(11-20) 


+ ksin\{i(if-+-^+2nCjg)-(^) 


k ^2 ^ 12 ^ 66 - 


'22 


A3(a) = iksin^acosa{-|^^-^(^+-^+2nCgg)(fs2+^2j 


2 Q 

+ C 2nC^^) }+ik cos\{-Cj^2 ^ 
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Equations (11-18) with (11-19,20) and the phase shifts a and 8 
given by (11-16) constitute the final form of the general solutions of the 
difference equations (11-13) . 

In Eqs. (11-19,20) we notice that when k 0 we have X^( 8 ) = X 2 ( 6 ) = 

- Y^(a) = 0. Namely the propagation of the normal stress (with phase 
shift 3 ) and the propagation of the shear stress (with phase shift a) are 
completely decoupled- This occurs when the waves are propagating only throiagh the 
thickness direction [27 ] . 

Impact Boundary Condition 

Boundary conditions for an impact can be described by any two conditions 

among u , v , a , and t and another two conditions from u»., v„, o„, 

® o’ o o’ o N N N and 

Tjj. For our present problem we have chosen a line impact by a normal stress 

along the axis (Figure 1), i- e., 

p 

o — ( 1 -cos -^^) ( 1 +cos^) : lx I < a and 0 < t < t 

o 4 t a — — — o 

o 

= 0: |x[>a or t>0 or t>t 

o 


T 

O 


'N 



(n-21) 


Hence, the boundary conditions for the present impact problem lead to "the 
following equation 





2 2 

V = {l+X 2 (e)Y 2 (a) }sin2aN.sin2$N + 2X2C6)Y2(ct) (,cos2aN-cos2eN-l) 

= X 2 (B)Y 2 (a) (cos2aN'cos26N-l) + sin2aN-sin28N 

= i{cos2BN.sin2aN-X2(3)Y2Ca)cos2aN-sin2aN} (11-23) 

= iX2(3){X2(S)Y2(a)sin2gN-cos2aN-cos26N-sin2aN} 

= X2(6){X2(S)Y2(a)sin2aN.sin2BN+cos23N-cos2oN-l} • 


Substituting the E|s into the general solution (ll-l8) we can find 

the complete solutions which satisfy the impact boundary conditions given 

by (11-21). In other words, for given values of k and s we first find 

the phase shift a and 3 from (11-15,16) and with these we can find 

solutions in integral transform from equations (11-18,23) which are the 

final solutions under impact. After U , V , T and I are calculated 

n’ n’ n n 

* Allowing the determinant of the coefficient matrix to vanish leads 

to dispersion relations of an K-layer plate, namely V{afi) = 0. Then, 
a and 3 are obtained from (Il-15,l6) which gives the complete 
dispersion relationships. 
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with a given impact function q, they can be inverted easily by means of the 
fast Fourier transform routine [3,20] to give the complete displacement and 
the stress fields after Impact. 

Tangential Normal Stress 

As discussed following Eq. (II-7),the tangential normal stress does not 
appear explicitly in the approximate equations of motion. Therefore, this 
component of the stress has to be calculated from the constitutive equation. 
Namely, 


'11 




'11 


n 


12 

c-n “ 1 + (v -V .) 1 < n < N 

11 n,l 2b n n-1 — — 


or after normalization and integral transform they are 


0^1 = -ikC^3_U^ + 


(11-24) 



= -ikC,.U + 
11 n 


N 


■‘=12‘VVl> 


1 < n < N . 


Then once the displacement field is computed the tangential normal stress 
can be obtained from the above equation and inverted. 



4. Some Numerical Results 


The analysis discussed in the previous section includes the transient 

propagation in all directions but suitable choices of impact time, impact 
radius, sizes of time and distance steps are essential to make good use 
of the fast Fourier transform. For example, if we take a large time 
increment with a relatively thin plate propagation through the thick- 
ness will not be seen. For this matter we have examined several different 
cases. 

Case 1: Longitudinal propagation 

Propagation of impact generated waves along the longitudinal direction 
is examined for an isotropic plate (steel plate: X = y = 1.2x10^ psi) 

employing a two-layer model. For these calculations we used an impact 
time t^ = 10 ysec, plate thickness A = 1 cm, and impact radius a = 4 cm. 
Some of the results at a few different time sequences are shown in Fig. 

6 a-f. 

In these figures we can see two distinct states of propagation and 
corresponding wave fronts : one for horizontal displacements (u) and 

longitudinal normal stresses another for vertical displacements 

(v) and shear stresses (t) . In other words, the initial signals of the 
horizontal displacements and longitudinal normal stresses propagate 
through the plate with longitudinal wave speed at amplitudes that are 
relatively small. But the major parts of their signals are due to a 
bending wave propagating with shear stresses and vertical displacements 
with a lower velocity. When the group velocities of these waves are 
calculated from the numerical results, they are about 5 mm/ysec and 
3 mm/ysec, respectively, while the phase velocities of the unbounded 
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medium of this material are = / (A+2ji) / p = 5.61 mm/ysec and 
Cg = /y/p = 3,25 mm/ysec. 

Case 2 : Propagation Through Thickness 

To examine the propagation through the thickness it is necessary to have 
a. sufficient number of layers in a plate. It is also essential to 
make the time step relatively small compared to the layer thickness. To 
do this we increase the thickness of the plate and the number of layers and 
reduce the Impact time. 

In Figs . 7 and 8 propagation of the transverse normal stress in the 
same plate (A = 4 cm, t^ = 2ysec, a = 40 cm; 4-layer model) is shown 
at different time sequences. As seen in Fig. 7»the transverse normal 
stress is initially compressive due to the impact and a compression wave 
propagates through the thickness. But later it becomes a tension wave 
after reflection from the free surface and the tension wave propagates 
back to the impact surface. In Fig. 8 we see the change of the transverse 
normal stress and the interlaminar shear stress with time for the same 
impact conditions as in Fig. 7. 

Similar results are also shown for the case of an anisotropic plate 
in Fig. 9 (55% graphite fibers-epoxy matrix, layup angle = 15®; A = 1 cm, 

= 2ysec, a = 2 cm; 8-layer model). Here we again notice a clear delay 
in time for waves to travel from one layer to the next one. Another important 
point is that the shape of the impact stress is relatively well preserved 
during the initial stage of propagation but changes immediately afterwards. 

The distortion of the shape becomes more serious with further propagation 
due to reflection, thus, showing the highly dispersive nature of the harmonic 

waves in the approximate plate theory. 
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When the group velocities are calculated from these resiolts, we find 6.32 
nim/ysec for the dilatation wave and 3-33 mm/psec for the shear wave in the 
case of the isotropic plate and 2.5 mm/ysec for the quasi-dilatation wave of 
the anisotropic plate. Their expected values are, respectively, 5. 61, 

3.25, and 2.36 mm/ysec. In other words, waves going through the thickness 
are traveling faster than expected. 


Case 3. Wave Surfaces 

In the previous two cases we examined the transient waves propagating 
dominantly along either the x^^ axis or through the thickness direction 
by suitable choices of the plate geometry and impact condition, 
now examine the combined effect, simultaneous propagation in both direc- 
tions. This effect is shown in Fig. 10 (isotropic plate; A = A cm, 
t^ =x Aysec, a = A cm; A-layer model) where the transverse normal stress 
generated from the line source due to impact not only spreads out in all 
directions but also reflects from the free surface. 

When the plate is anisotropic, the situation is more complex in the 
sense that waves are neither dilatation nor shear but they are coupled 
together (now called quasi-dilational or quasi-shear waves) . Due to the 
coupling, phase velocities of the anisotropic wave vary from one direction 
to another, resulting in complicated shapes for the velocity surfaces and 
wave fronts [2], For an ansiotropic plate (made of 55% graphite fiber- 
epoxy matrix with la3nip angle A5**) the velocity surfaces and the wave 
surfaces are shown in Fig. 11. The stress state at 10 ysec after the 
impact on the same plate (A = A cm, t^ = Aysec, a = 2 cm; 8-layer model) 
with the corresponding wave fronts are shown in Fig. 12a. In the 
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propagation of the quasi-longitudinal wave we notice that the longi- 
tudinal propagation is well bounded by the quasi-dilatational wave surface 
but the transverse propagation is not. The shear wave is not bounded by 
the quasi-shear wave front in either direction. 

This interesting phenomenon of higher propagation speeds through the 
thickness is related to the dispersion relationship at short wave length 
limits; it is discussed in the next section. 



5. Correction Factor and Conclusion 


Correction Factor 

According to the present model of a multilayer plate, one of the 
antisymmetric modes of the dispersion relation^ips approaches the shear speed 
when the wave length becomes shorter and shorter, as mentioned in Section 
2. It is well understood that such a limit is incorrect, i.e., in the 
limit of short wave length there should be a Rayleigh wave instead of a 
shear wave. Such a discrepancy can be eliminated by introduction of 
proper correction factors, as shown by Mindlin and Medick P-8 ] . Correction 
factors can be found by examining either the large wave number limit or 
the cut-off frequencies of both the exact theory and the present approximate 
theory. Since these two ways lead us to the same results we will 
find the factors by matching the cut-off frequencies of the two theories. 

The cut-off frequencies of the exact theory for an isotropic plate 
can be obtained from the well-knora Rayleigh-Lam’s equation. The 
lowest cut-off frequencies are /(A+2p)/p for the symmetric mode and 
/ y /p for the antisymmetric mode. The corresponding cut-off frequencies 
of our approximate theory are — and — / 3cgg/p .. Hence, we notice 

that replacing 0^2 ^22^ ^66 ^66^ makes the two 

theories have the same two lowest cut-off frequencies. Furthermore the 
shear wave observed in the short wave length limit of the present approx- 
imation becomes a wave with a speed of — ^ /y/p, i.e., the Rayleigh wave. 

/12 

Another important consequence of the correction factor is to reduce 
propagation speeds through the thickness, which are related to '^^2lJ^ 

* 

The lowest two cut-off frequencies are found from Eq. (II-ll) and they 
are independent of the layer number in the plate under investigation. 
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and /c- ,/p , with a factor of ir/ »^T2 . Propagation of the maximum value 
bo 

of the interlaminar normal stress through the thickness is examined with 
and without correction factors and the results are shown in Fig. 13. 

Without the correction factor the propagation speed in a composite plate 
is roughly about 2.60 mm/ysec obtained from the numerical results used in Fig. 12. 

When the same plate is subjected to identical impact conditions this 
reduces to about 2.41 mm/ysec with the correction factor. Comparing this 
with the group velocity in an unbounded composite space (= 2.36 mm/ysec) 
the agreement of the present approximate theory is remarkable. Similar 
results are also observed in the case of shear and quasi-shear waves. When 
these correction factors are Introduced in the previous cases, shown in Figs. 

8 , 9, and 12j all the signals propagating through the thickness are now 
well bounded within the corresponding wave fronts, as shown in Fig. 12b 
and from this we can notice the importance of the correction factors. 

Discussion and Computation Time 

It is interesting to compare the computation time of this model with 
some other methods, such as the finite element method or the finite difference 
method. In the case of an 8-layer anisotropic plate model, from which Figs. 9 
and 12 are produced, we have 

9 steps along the thickness; 8-layer model; 

32 step along the Xj^ direction; 64 points are used in pratice but 

only half of them are useful because 
of the symmetry of the problem, 

32 steps in time; 

2 displacement components at each point. 

Therefore the total number of the unknowns, which are the basic unknowns 
either in case of the finite difference or finite element methods, is 



18,432. After these primary variables are calculated, 27,648 secondary 
variables (three stress components at each points) have to be calculated 
again. According to our present model all these processes require only 
200 K of computer space without using magnetic tapes or any kind of 
additional storage space and only 1 minute 6 seconds for CPU time in the IBM 
370-168 model including compiling, linkage editing, 1/0 and execution. 

Conclusion 

The present theory is a generalization of Mindlin's approximate plate 
theory applied to a multilayer plate under an Impact. By combined use of the 
finite difference technique in the thickness direction and the fast Fourier 
transform in the plane of plate and time, this model can be very useful 
for the study of wave propagation in a composite plate under impact forces. 
However, reasonable attention in usage of the fast Fourier transform is 
required to avoid spurious data. From the limited numerical data obtained 
from this model it appears that the anisotropy in the plate will lead to a 
considerable interlaminar shear which might lead to ply bonding failures. 

The model also shows that for short enough impact times, an interlaminar 
tension can develop as one would expect, which might also account for 
interlaminar ply failure. 



III. IMPACT OF A COMPOSITE PLATE WITH AN INTERLAMINAR DAMPING LAYER 


1. Description of Problem 
Geometry of Plate 

As an extension of the multilayer plate discussed in Chapter II 

ve nov examine the impact and the consequent stress 

wave propagation in a composite plate with viscoelastic damping layers. 
Possible models for damping mechanisms in plates are shown in Fig. 14. 

We will formulate a model made of an alternating 

number of elastic and viscoelastic layers, as shown in Fig. l4-c . 

As long as the layer structure of the plate is periodic, the main part of 
the analysis in Chapter II for an elastic plate is valid with additional 
equations for viscoelastic layers. 

Viscoelastic Property of Elastomer 

The mechanical properties of an elastomer are usually expressed in 
terms of a complex modulus depending on the frequency, i.e., 

G*(aj) = G'(o)) + iG"(u>) . (III-l) 

With this the constitutive equation is written as 


(oj) = G (u))e^j(o)) 


(III-2) 


in the frequency space where ^ are respectively the 

A 

Fourier transforms of a^. and e.. in time [29]. 

ij ij 


For (III-2) the Fourier transform is defined as 


f (to) = f f(t)e'^“^dt . 


-a 


TrTr«^ATj PAGE IS 
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The constitutive relation (III-2) with the complex modulus (III-'l) 
implies the following constitutive .equation in a time space: 

t 

o(x,t) = / G(t-T)e(x,T)dT (III-3) 

where the relaxation function G(t) is related with the complex modulus 
•k 

G (u)) as 

00 

G*(m) = f G(t)e"^“*^dt . (III-4) 

o 

k 

Therefore, when the complex modulus G (to) is obtained by experiments, 
usually by means of harmonic excitation of strain, the relaxation function 
G(t) can be found by inversion of equation (III-4) . 

The viscoelastic property of the elastomer under consideration has 
been extensively investigated (s.g. [19]) its 

complex modulus is given in Fig. 15. This complex modulus can be reasonably 
well described by a three parameter equation as 

G'(to) = a - . (HI-5) 

(0 +c 

These three parameters are obtained from another set of parameters: the 

maximum values of G'(o)) when or*®, the maximum value of G" (o>) /G ' (oj) 

and the at which G" (to) /G ' (oj) becomes the maximum. Therefore »if 

we characterize the complex modulus by proper choice of G'(oj), to and 

o 

the maximum of G"(o) )/G'(oj ), the relaxation functions are completely 

o o 


described. 



2. Formulation 


Elastic Layer 

In Fig. 16 a typical viscoelastic layer (nth) is shown between two 
adjacent elastic layers (nth and (n+l)th) with appropriate discretiza- 
tion. The approximate equations of motion for the nth elastic layer 
given by Eq. (II-7) are still valid. But remembering the new discretizing 
notation in Fig. 16 we now have to replace ( and ^ 1 ^ 

and ( )"*" 1 5 respectively. The results are 
n-i 


p (u +u"^ ^ , (u +u^ , ) ^ , H — ^(v -v^ . ) . +'r(f -'t'*' , ) 

n n-1 11 n u-1 ,11 b n n-1 ,1 b n n-1 




12 . — + » 22 f + \ \ \ ~ * \ 

— f— (u+u ,) . T'^v^-v^ i)+v:(o+o„ i) 1 

b n n-i ,i n n-1 b n n-1 n n-1 ,1 


3 , + X . , - + 


P 


3c 3c 

- / - + X 66, - + X 66, -X. + X 

‘=ll<VVl> .11 --2-'%-“o-l) <''n-^n-l> ,1 


(III-6) 


. 12,— + X .3,-.+ X 
+ — ^(o -a -) . +■'■ i) 
22 n n-1 ,1 b n n-1 


P(v~+v^ x) 
n n-1 


‘^66^i^vCi^ ,1-"^ Vvi^ .11^ b^ vvP 


Viscoelastic Layer 

Since the thickness of the elastomer is thin compared with the elastic 

layer, we can assume that the stress field is uniform through the thickness 

- + — + 

of the elastomer. In other words, we have a = a = a and t = t = t 

n n n n n n 

for (III-6) . Therefore, the following compatibility conditions for the 
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elastomer can be obtained immediately: 


'12 (n) - 4 + 2D<V“n' 

(III-7) 

'22(n) - 

where D is the thickness of the elastomer. 

We further assume that the dissipation is mostly due to shear motion, 
( i.e., that the normal component of the continuous traction vector is 

transmitted through the viscoelastic layer purely elastically. There- 
fore, by combining (III-7) with (III-3) we find 


a 

n 


E. + 

— (v -V ) 
D n n 


T 

n 


/ G(t-T) •^(v‘!‘+v“) +'^(%-u j }dT 
2 ox, n n D n n 
_« 1 


(III-8) 


These two equations and four more from Eq. (III-6) are the complete 

equations needed to solve the impact on a composite plate with elastomer. 

For a plate made of N elastic layers and (N-1) viscoelastic layers 

Eq. (III-6) provides 4N equations and (III-8) gives 2 (N-1) equations. 

Since the total number of the unknowns are now 6N+2 (u , v , a , x ; 

o o o o 

-- + + -- + + 


Ui* v^, v^, o^, T^; 

u^, v^, a^, Tjj) we can solve this system of equations with four additional 
conditions supplied by the suitable boundary conditions. 


Here we notice that the governing equations are now a set of six 
difference-integro-partial differential equations. These equations can be 
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reduced to difference equations after appropriate integral transforms and 
the resulting difference equations can be handled rather simply, as in the 
previous chapter. 

3. Numerical Results and Discussion 
Impact on Plate 

For the report ve examine the impact on a plate consisting 
of two elastic layers and a viscoelastic layer»as shown in Fig. 17>with 
an impact function 

2 

a = P {l-( — ) } sin — ; lx, I < a and 0 < t < t 

o o a t 1 o 

o 

(III-9) 

= 0 ; |xj^| > a or t > t^ , or t < 0 

with all other stress components vanishing on both surfaces of the plate. 

Now by putting n = 1 and 2 into Eq. (III-6) we have eight equations and 
two more equations are obtained from Eq. (III-8) . We again normalize 
these equations and take the integral transform, as in Chapter II. The 


resulting equations are : 
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-(fs^+^^k^)(u7+U )-CT,ik(v7-V )+T = 0 

a 1 o J-iC 1 o 1 

-<£=^+|cj^lc2)(v-4.5^)-Cjjik(§^-5_,)4j = 


-(fs^+^T,k^+^3C,,)(u"-U )+3C ik(v“+V ) - jikf +3T 

' A 11 b 66 1 o 66 1 o 0^2 All 


C 22 A 0 


2 n 2 - i ;i-i. ^ 

c ^ ^ 

-(fs^ +^U‘'^ *f^‘=66> * C^ ^VA 


h ■ W^) 


§2 - g(sH^(0^4-)-f(5;;4-)} 


where G(s) is the Laplace transfora of the relaxation function G(t) 
with respect to x = t/T^ and we have used the boundary conditions 

= T 2 = ® t^ke above 10 equations we can find 10 unknowns 

A A I A I A A A A A 

(U , V , U , V , u. , V,, I, , T- ; U-, V ) with given impact function Z 
OOlXiiXlZ/ o 

and once these are calculated the displacement and the stress fields can 
be computed by inversions of the integral transforms by means of the FFT 
alogorithm. 
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Numerical Results 

For the present computation we have used the Young's modulus 

12 

A 42 41*10 

E = ,7*10 psi and the shear modulus G'((d) = .817*19 — r — ^ for 

3*10 -ho 

the elastomer where oo is given in hertz. The G'(o)) in this case implies 

that G'(“) = .817*10^ psi and max(G"(to)/G' (to)) =3.3 at u = 800 Hz. 

o 

The propagation of stress wave in this case is quite similar to that 
of the composite plate without an elastomer layer except the peak values 
of the interlaminar stress. Values of the peak stress with different 
thickness of the elastomer layer are plotted with those of the purely elas- 
tic plate in Fig. l8. As we can see in this figure the interlaminar shear 
stress has increased by a small amount while a reduction of the normal 
stress is considerable when the elastomer layer becomes thicker and thicker. 
From this result it is obvious that the reduction of the normal component 
of stress can be achieved by introducing such a soft and energy-dissipating 
elastomer layer. 

Discussion 

In addition to the simple reduction of the normal stress it is also 

observed that the amount of reduction increases with the value of G"(iu)/G' (m) 

and the location of u at which G"(u))/G' (u)) becomes the maximum value. 

o 

In other words, we can make the dissipation effect more serious by choosing 

an elastomer whose G"(aj)/G' (m) becomes maximum at m around which the 

o 

most of the impact energy is carried out. 

It is also believed that a further dissipation effect will be possible 
if we make the transmission of the normal stress viscoelastic across the 
elastomer layer, which we have assumed is elastic for this report. 



IV. IMPACT ON A PLATE WITH A CRACK 


1. Introduction 

When the impact stress is low, the impact is elastic and the 
stresses in the plate can be described by elastic wave propagation. When 
the stress is increased beyond a certain limit then the impact damage 
occurs. Elastic-plastic impact is complicated for two reasons, namely, 
unloading and loading must be treated differently, and the strain rate effect 
[30] must be included. If the impact stress is increased further to a 
certain level where the induced stress is higher than the strength of a 
target material then penetration begins to occur. In this limit the 
target material sometimes behaves as a fluid and such a state of impact is ' 
known as a hydraulic impact [31]. Another failure mode is the occurrence of 
interlaminar cracks. 

Investigation of the stress state in solids with cracks falls in the 
category of so-called fracture mechanics and has been under an extensive 
scrutiny since the famous enunciation by Griffith [32]. Presence of cracks 
inside a material usually leads us to a mixed boundary value problem and 
only a limited class of problems can be solved [33,34] . In the case of 
dynamic loading the problem becomes more difficult due to the scattering 
of the stress wave by the crack [21-24]. In this report we will 
formulate the problem of a plate with a crack which is subject to a 
d 3 mamic loading. 

Our original goal was to study the effect of interlaminar cracks in 
composite plates in response to impact loads. Debugging problems in 
other parts of this report, however, used valuable time originally set aside 
for this problem. The following section is an attempt to illustrate the 
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use of the Mindlin plate theory for the study of interlaminar cracks and 
to point out the mathematical difficulties that must be overcome in solving 


the problem. 



2. Formulation 


Description of Problem 

The plate under consideration has a crack on the midplane running 
from x^= -h to +h as shown in Fig. 19. Stress can be applied either 
on the surface of the plate or on the crack surface. In the former 
case the crack surfaces can be in contact and the boundary conditions 
become more complex due to the partial continuity of stresses and dis- 
placements during the contact. For the present report to illustrate the 
mathematical difficulties we assume that the crack surface is subject to 
a known compressive impact. 


Governing Equation and Boundary Conditions 

We can formulate this crack problem by assuming that the lower and 
the upper half plates are made of a number of layers but for 
simplicity we consider the plate to consist of two identical layers and the 
crack to be present on the interface of these two layers. Following the 
notation shown in Fig. 19 we have the governing equations identical to 
Eq. (III-6) with n = 1 and 2. The boundary condition requires that 
both plate surfaces remain traction free. The crack surface is subject 
to a prescribed impact condition while the displacement and stress are 
continuous along the layer boundary outside the crack. Namely, we have 


= T 


02 = t2 = 0 


'1 

+ 


= o. 


-P (x ,t) 
o i 


= T, = 0 


r |x| < h 


u, = u. 


= O, 


+ — 
Vi = v^ 


X > h 


(IV-1) 
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Due to the twofold symmetry of the problem we now have = U 2 > 

+ — + — + ~ 

= “^ 2 , "'’i ~ ® '^l ~ '^1 ~ '^l’ ~"'^1 ~ ''^1 ~ ^ 1 ’ 

Thus^ the eight equations obtained from Eq. (III-6) are now reduced to 




p(v -V ) 
1 o 


p(u^-ii ) 
± o 


^‘l2, . , ^' 22 , > . 3(J 

D 


(IV-2) 


= c,,(u,-u_) 


3c 3c c 

-(u -u^) -— j— (v,+v„) T +■ 


--I *1 y -1 -I ^ V“-i *^ / 1_ W ,, • V / - • w T 

11 1 o ,11 ^2 ' 1 o b 1 o ,1 C 22 >1 


pCvi-f^o^ = ‘'66^i^V“o\l‘"^V^\ll^'"f 


and the boundary condition is now 


a = -P (x, , t^ X < h 

o 1 I 


V, = 0 


lx| > h 


along X 2 = 0 


(IV-3) 


Dual Integral Equation 

We now normalize the governing equation (IV-2) and take the integral 
transform. Then we have 


(fs^+^llk^), ^12^^ 

0 

9 

0 

-3Ci2ik , (fs +"b^22^’ 

0 

9 

0 

0 

0 

"66""^ 

9 
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0 
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and these can be solved for U , U, , V , and V, in terms of E. Since 

o 1 o i 

the mixed boundary conditions are given by a and v^ we solve as 


= K(s,k)E 


(IV-5) 


with 


K(s,k) 


lr3,£ 2 b_ ,2^ .If 


12 


22 




A = Det 


(fs^+|Cj^^k^) , C^^IV. 


-3C3_2ik 


’ (fs +-^-^22^ 


(IV-6) 


B 


Det 




/•fg2 1^2 3A_ V 

^^11^ T^66^ 


te 2 ,2. 

+^66^ > 

-3C,,ik 

66 


Next we take the Inverse transform of E and V^, and apply the mixed 
boundary condition given in Eq. (IV-3). Since the boundary conditions are 
for all times t > 0 we only take the inverse Fourier transform to apply 
the boundary conditions, i.e.. 


00 

E(n,s) = — / E e“^‘^’^dk 



00 

V (n,s) = — / K(s,k) E e'^^’^dk 
-00 


(IV-7) 


Application of the boundary condition given by Eq. (IV-3) results in the 
following integral equation: 
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-p^(n,s) 



^ -ikn , , 
Z e dk 


n| < h/A 


00 

0 = — / K(s,k) Z e~^^’^dk : |n| > h/A 

— oo 


(IV-8) 


for an unknown function E. 

3. Discussion 

The integral equations of the type given in Eq. (IV-8) are known as 
dual integral equations, each of which has its own region of application 
and occur in mixed boundary value problems [35]. There are a number of 
ways to solve this type of integral equations, such as by reduction to a 
single Fredholm integral equation or hy using the Wiener-Hopf technique [36]. 
Finding the solution depends on the kernel and in general it is rather 
difficult to do except for some speciaJ. cases such as for Bessel kernels or 
trigonometric kernels. 

Once the unknown function Z is determined the other variables 

can be computed by solving the algebraic equation (IV-4) 
and the complete displacement can be found by inversions of the integral 
transform. 

The problem formulated in this chapter is the simplest impact 
problem in that the contact of the crack surface does not occur and that 
it has a twofold symmetry. But it is expected that the critical response 
of the plate, particularly the stress field near the crack, can be a guide- 
line for a more complex problem. 



V. CONCLUSION AND RECOMMENDED RESEARCH 


The present theory is a generalization of Mindlin's approximate 
theory of plate applied to a multilayer plate under impact. By combined 
use of the finite difference technique in the thickness direction 
and integral transforms this model has been shown to be very effective 
for wave propagation analyses . 

This model is extended to examine the effects of an elastomer layer 
between elastic layers of the plate. The reduction of interlaminar normal 
stress is signficant due to the damping layer but further investigation 
seems necessary to determine the nature of the reduction. 

The presence of a crack in the plate has been form^ilated. The re- 
sulting equations are given by dual integral equations which, as in many 
cases, are rather diffictilt to solve. 

The basic idea of the periodic structiare of the multilayer plate, 
where the governing equations are derived for each layer and given by a 
set of difference-differential equations, may be useful to handle dif- 
ferent types of problems, such as heat conduction and thermoelastic 
problems in composite plates . 


- 1+3 - 




References 


[l] Mindlin, R.D., - "High Frequency Vibrations of Crystal Plates", 

Quart. Anpl. Math ., Vol. 19 , p. 51-61 (1961). 

[2] Moon, F.C., - "Wave Surfaces Due to Impact on Anisotropic Plates", 

J. Composite Materials , Vol. 6 , p. 62 (1972). 

[ 3 ] Moon, F.C. and Kang, C.K., - "Analysis of Edge Impact Stresses in 
Composite Plates", TR to NASA Lewis Research Center, Princeton 
University (1974) . 

[li] Moon, F.C., - "One Dimensional Transient Waves in Anisotropic Plate", 

J. Appl. Mech. , Trans . ASME , p. 485 (1973). 

[ 5 ] Moon, F.C., Kim, B.S. and Fang - Landau, S.R., - "Impact of Composite 
Plates: Analysis of Stresses and Forces", TR to NASA Lewis Research 

Center, Princeton University (1976). 

[g] White, J.E. and Angona, F.A., - "Elastic Wave Velocities in Laminated 
Media", J. Acoust. Soc. Am. , Vol. 27, p. 310 (1955). 

[ 7 ] Sun, C.T., Achenbach, J.D. and Herrmann, G. , - "Continuum Theory for 
a Laminated Media", J . Appl . Mech . . Trans , ASME , p. 467 (1968). 

[ 8 ] Nayfeh, A.H. and Gurtman, G.A., - "A Continuum Approach to the 
Propagation of Shear Waves in Laminated Wave Guides", J . Appl . Mech . . 
Trans. ASME , p. 106 (1974). 

[ 9 ] Hegemier, G.A., Gurtman, G.A. and Nayfeh, A.M. , - "A Continuum 
Mixture Theory of Wave Propagation in Laminated and Fiber Reinforced 
Composites", Int. J. Solids Struc. , Vol. 9, p. 395 (1973). 

[10] Chao, T. and Lee, P.C.Y., - "Discrete Continuum Theory for Periodically 
Layered Composite Materials", J. Acoust. Soc. Am., Vol. 57, p. 78 
(1975). 

[11] Lee, P.Y.C. and Nikodem, Z., - "An Approximate Theory for High Frequency 
Vibrations of Elastic Plates", Int. J. Solids Struc., Vol. 8 , p. 581 
(1972). 

[12] Dong, S.B. and Nelson, R.B., - "On Natural Vibrations and Waves in 
Laminated Orthotropic Plates", J. Appl. Mech., Trans. ASME, p. 739 
(1972). 

[ 23 ] Hegemier, G.A. and Nayfeh, A.H. , - "A Continuum Theory for Wave 

Propagation in Laminated Composites, Case 1: Propagation Normal to 

the Laminates", J . Appl . Mech . , Trans. ASME , p. 503 (1973). 

Hegemier, G.A. and Nayfeh, A.H., - "A Continuum Theory for Wave 
Propagation in Laminated Composites, Case 2: Propagation Parallel 

to the Laminates" J. Elasticity . Vol. 3, p. 125-140 (1973). 


~ kk - 



[15] Cooley, J.W., Lewis, P.A.W. and Welch, P.D., - "The Fast Fourier 
Transform Algorithm: Programming Considerations in the Calculation 
of since, cosine and Laplace Transforms", J. Sound Vib . , Vol. 12, 
p. 315-337 (1970). 

[16] Brillouin, L. - Wave Propagation in Periodic Structures , Dover, N.Y. (19U6). 

[17] Thomson, W.T., - Vibration Theory and Application , Prentice Hall, 

N.J. (1965). 

[18] Mindlin, R.D. and Medlck, M.A., - "Extensional Vibrations of Elastic 
Plates", J. Appl. Mech. , Trans. ASME , p. 561 (1959). 

[19] Yan, M.J. and Dowell, E.H., - "High Damping Measurements and 
Preliminary Evaluation of an Equation for Constrained Layer Damping", 

AlAA Journal, Vol. 11, p. 388 (1973). 

[20] Nakra, B.C. and Grootenhuis, P., - "Structural Damping Using a 
Four Layer Sandwich", J. Eng. Ind. , Trans. ASME , p. 81 (1972). 

[21] Achenback, J.D., - "Crack Propagation Generated By a Horizontally 
Polarized Shear Wave", J. Mech. Phys. Solids , Vol. 18, p. 245 (1970). 

[22] Achenback, J.D. and Nuismer, R. , - "Fracture Generated by a Dilatational 
Wave", Int. J. Frac. Mech. , Vol. 7, p. 77 (1971). 

[23] Nuismer, R.J. and Achenback, J.D., - "D3nnamically Induced Fracture", 

J. Mech. Phys. Solids , Vol. 20, p. 203 (1972). 

[2U] Keer, L.M. and Luong, W.C., - "Diffraction of Waves and Stress Intensity 
Factors in a Cracked Layered Composite", J. Acoust. Soc. Am., Vol. 56, 
p. 1681 (1974). 

[25] Levy, H. and Lessman, F., - Finite Difference Equations , MacMillan Co. 
(1961). 

[26] Kim, B.S. and Moon, F., - "Impact on a Laminated Half Space", to be 
published . 

[27] Premont, E.J. and Stubenrauch, K.R., - "Impact Resistance of Com- 
posite Fan Blades", Tech. Rep. to NASA Lewis Research Center, NASA 
CR-134515, Pratt and Whitney Aircraft (1973). 

[28] Novak, R.C., - "Multi-Fiber Composites", Tech. Rep. to NASA Lewis 
Research Center, NASA OR-135062, United Technology (1976). 

[29] Christensen, R.M. , - Theory of Viscoelasticity , Academic Press, N.Y. 

(1971). 

[30] Christescu, N., - Dynamic Plasticity , North-Holland , New York (1967). 

[31] Birkhoff, A. ^ - "Explosives with Lined Cavities", J. Appl. Phys . 

Vol. 19, p. 563 (1948). 



- k6 - 


[ 32 I Griffith, A., - "Phenomenon of Rupture and Flaw in the Solids", 

Phil. Trans. Roy. Soc. (London), Vol- A221, p. 163 (1921). 

[ 33 ] hiebowitz, H. , (Ed.) - Fracture , Vol. 2, Academic Press, N.Y. (1968). 

Sneddon, I.N. and Lowengrub, M. - Crack Problems in the Classical 
Theory of Elasticity , John Wiley and Sons, New York (1969). 

[ 35 ] Sneddon, I.N., - Mixed Boundary Value Problems in Potential Theory , 
North-Holland, Amsterdam (1966). 

[36] Noble, B., - Method Based On the Wiener-Hopf Technique , Pergamon 
Press, New York (1958). 



Fisrures 


I. Composite Plate and Layer 

2 a,b. Dispersion Relationship and Phase Velocity for Isotropic Plate: 

One-layer Model (X = y) . 

3 a,b. Dispersion Relationship and Phase Velocity for Composite Plate: 

One-layer Model (55% Graphite Fiber-Epoxy Matrix, Layup Angle 45“). 

4 a,b. Dispersion Relationship and Phase Velocity for Isotropic Plate: 

Two-layer Model (X = y) . 

5 a,b. Dispersion Relationship and Phase Velocity for Composite Plate: 

Two-layer Model (55% Graphite Fiber-Epoxy Matrix, Layup Angle 45®). 

6 a~f. Longitudinal Propagation Impact Stress in Isotropic Plate: 

Two-layer Model (X = y = 1.2*107 psi; A = 1 cm, t =10 ysec, 
a = 4 cm) . ° 

7. Transverse Propagation of Normal Stress in Isotropic Plate: 

4-layer Model (X = y = 1.2*107 psi; A = 4 cm, t = 2 ysec, a = 40 cm). 

o 

8. Transverse Propagation of Normal and Shear Stress in Isotropic Plate: 
(Same as in Fig. 7). 

9. Transverse Propagation of Normal Stress in Composite Plate: 

8-layer Model (55% Graphite Epoxy-Fiber Matrix, Layup Angle 15®; 

A = 1 cm, t = 2 ysec, a = 2 cm), 
o 

10. Two Dimensional Propagation of Normal Stress in Isotropic Plate: 

4-layer Model (X = y = 1.2*10^ psi; A = 4 cm, t^ = 4 ysec, a = 4 cm). 

II. Velocity Surface and Wave Surface of Composite Plate (55% Graphite 
Fiber-Epoxy Matrix, Layup Angle 45®). 

12a, b. Comparison of Theoretical Wave Front and Numerical Wave Front of 
Composite Plate 10 ysec after Impact: 55% Graphite Fiber-Epoxy 

Matrix (For Numerical Results; 8-layer Model, A = 4 cm, t =4 ysec, 
a = 2 cm). Without and with Correction Factors. 

13. Effect of Correction Factor or Transverse Propagation of Peak 

Value of Normal Stress : 8-layer Model (55% Graphite Fiber-Epoxy 

Matrix, Layup Angle 45®; A = 4 cm, t^ = 4 ysec, a = 2 cm). 

14. Viscoelastic Impact Energy Absorbing Models. 

15a, b. Complex Modulus of Elastomer. 

16. Plate with Viscoelastic Layers. 

17. Impact of Plate Made of 2 Elastic Layers and a Viscoelastic Layer. 
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18. Peak Value of Interlaminar Stress Vs Elastomer thickness (Two 

Elastic Layers and a Viscoelastic Layer: (A = 1 cm, t =10 ysec, 
a = 4 cm) . ° 

19. Composite Plate with Crack. 
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Fi^e 2. Dispersion Relationship and Phase Velocity of Isotropic Plate 
One-Layer Model (A = p , Poisson's Ration = 1/4) ' 
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Figure 7 . Transverse Propagation of Normal Stress in Isotropic 
Plate: 4-layer Model (X- = y = 1.2*107 psi; 

A = 4 cm, t = 2 ysec, a = 40 cm) 











Impact surface 


I \/ C— 2.5mmy^sec 


V 


Free surface 


Figure 9 Transverse propagation of normal stress in 

composite plate; 8-layer Model (55% Graphite 
Fiber - Epoxy Matrix, ±15° Layup; A = 1 cm, 
t_ = 2 usec, a = 2cm) 
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Figure 11, Velocity Surface and Wave Surface of 
Composite Plate C557o Graphite Fiber- 
Epoxy Matrix, Layup Angle 45*) 
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Figure 12a .Wave front 1C nsec after impact (without correction factor) 
(557o Graphite Fiber-Epoxy Matrix, ±45° Layup; A = 4 cm, 

8- layer Model; t^ = 4 ysec, a = 2 cm) 
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Figure 12b.Wave front 10 fisec after impact (with correction factor) 
(55% Graphite Fiber-Epoxy Matrix, ±45° Layup; A = 4 cm, 
8- layer Model; t = 4 ysec, a = 2 cm) 
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FIGURE 15a SHEAR MODULUS AND LOSS TANGENT AT T ®F vs 
FREQUENCY FOR DUPONT LR3-604 
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Figure 19. Composite Plate with Crack 



AFPEITDIX A FLOW CHART 

In this flow chart and program, U(I,J), V(I,J), TAU(I,J) 
SIGMA (I,J) and SIGMA1(I,J) represent U, V, T and Z in 
Eq. (11-17,18) and integral transform of 


(i) Define the variables (U,V,TAU, SIGMA) for displacement and 
stress fields in matrix form 

(ii) Define other working matrices 

(iii) Supply the core space for the matrices 

i ^ 

(i) Read the elastic moduli of the composite plate (in PSI) 
and write them 

(ii) Normalize them by cgg 

(iii) Introduce correction factors for C 22 Cgg 

(iv) Calculate Egg (^£6 unit) and store 


1 


(i) 

Read and 

write the data for geometry of composite plate. 


NLAYER 

: Number of layers in plate 


DELTA: 

Thickness of 

the plate (in cm) 


RHO: 

3 

Density of the composite (in gr/cm ) 

(ii) 

Calculate 

the B (half of 

the layer thickness) 

(iii) 

Convert all quantities in 

MKS unit 



(i) 

Read and 

write all the impact data 


NA: 

Impact radius 

by integer multiple of plate 



thickness (A 

= NA* DELTA) 


TAUO: 

Impact time (in second) 

(ii) 

Read and 

write all the data of integral steps 


NX: 

Total step number in space 


NT: 

Total step number in time 


NXIMP: 

Step number in impact radius 


NTIMP: 

Step number in impact time 

(iii) 

Calculate 

step size in time and space 


Space : 

DX = A/NXIMP 



Time: 

DT = TAUO/NTIMP 
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(i) Calculate normalization units 

Space: UNITX = DELTA (Thickness of the plate) 

Time: UNITT = A/ /Egg/RHO (Time refuired for the quasi- 

shear wave to travel the impact radius) 

(ii) Normalized all quantities by UNITT and UNITX 

(iii) Calculate integral limits (m^ and k^) in Fast Fourier Transform 


Calculate QO(I,J), CBX(I,J), CAX(I,J), XIBX(I,J), YIAX(I,J) 

Y3AX(I,J), over a half of the range of integration 
(I = 1 - NT, J = 1, NX/2). 

QO: Impact function given in Eq. (11-22) 

CBX, CAX: cosB and cosa in Eq. (11-16,17) by DPHASE 

XIBX, X2BX ...: Xi(g), X 2 (B) in Eq. (11-18, 19) by DELL 

YlAX, Y2AX ...: Yi(a), Y2M in Eq. (11-18,19) by DELL 


Invert (see Block A) and check the impact function with QO 





- - 



N = 0 (Impact Surface) 


Over the range of data 

generation (I=1-NT, J=l- NX/2) 

Calculate 

a 

and 8 from cosa and cos 8 

Calculate 


. V-^y .. in Eq. (11-23) from DET 

Calculate 

cos 2nS. cos 2na, sin 2nS and sin 2na 

Calculate 

U(I,J), V(I,J), TAU(I,J), and SIGMA(I,J) 
in transformed space by Eq. (11-18,23) 


Invert (see Block A) U , V , T and Z for stress and displacement 
u n’ n n n 

field 


Calculate the tangential normal stress SIGMA1(I,J) by Eq. (11-24) 
in integral transformed space and invert (see Block A) 


N= NLAYER?^-\^ 
(Have we finished calculation 
''^^'\for all layers?) 


N = N+1 

(Go to next layer and 
repeat calculations) 


STOP 









Block A Inversion 


(i) Data xx(I,J) in integral transformed space are generated 

for a half of the inverse transform range: I = 1 ~ NT, 

J = 1 -• NX/2 

(ii) Generate full data over J = 1 - NX by FLIP 

Symmetric flip: V(I,J), SIGMA(I,J), SIGMA1(I, J),QO(l,j ) 

Antisymmetric flip: U(I,J), TAU(I,J) 

(iii) Invert them for displacement and stress fields by FOURT 

(iv) Take care of the coordinate shifts and multiplication 
factors in FOURT by FACT 

(v) Print out by MAP 
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APPENDIX B SAMPLE COMPUTER DECK - QUALITy 
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3d S 

24 


Data Card 4 1 

1 

/“ 

c. 




1 

Data Card 3 | 

7 

c* 

*«• 

i. 


7. 44 

Data Card 2 ;j 

1 



PROGRAM IN SOURCE DECK (see Appendix C) 


IC FE3RTGCLG 


REGJ.DH=£07K j LIHES-SK ? CLhSS=C 


JOB CARD 
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I I ] 4 S t ; I 9 I3 n :: O U 19 It n II 19 73 n 29 2124 n » ZI 2I79U II U ]] >4 nu 97 19 19 40 41 42 4144 49 4t II 11)190 91 92 91 94 99 99 91 99 99 10 91 12 92 14 99 Cl II II El 10 n n II 14 19 It n II 19 10 

1 1 11 1 1 11 11 1 1 1 1 1 1 11 1 1 n 1 1 11 n 1 11 1 1 n n 1 1 1 1 11 1 1 n 1 1 11 1 1 n 1 1 1 n 1 11 1 1 1 11 1 11 n n 11 1 1 1 
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APPENDIX C LISTING OF PROGRAM AMD SAMPLE OUTPUT 


1 RELEASE 2.0 


MAIN 


DATE = 77139 


21/11/39 


5r s|t 4s ** Jl* * * 4= 4: ^ sit ^ ’S' ***=!!=<=* sf: sj: sf: 5!: 4: * sf: 5|: * =i= =!= ’l' *** ^ Jit 


THIS PROGRAM CALCULATES THE TRANSIENT PROPAGATION OF STRESS WAVE 
IN A LAMINATED COMPOSITE PLATE DUE TO A NORMAL IMPACT. 


4£4'4=4::ir: 


THE PRESENT PROGRAM IS A PART OF THE RESEARCH PROJECT OF 
PROFESSOR FRANCIS C. MOON 

DEPARTMENT OF THEORETICAL AND APPLIED MECHANICS 
CORNELL UNIVERSITY 

AND SUPPORTED BY NASA-LEWIS RESEARCH CENTER. 


4: 4= 4= 4: 4: 4t 4= 4« 4: 4s 4« 4t 4: 4: 4s 4= 4s 4: * 4: 4= 4= 4! 4: 4= 4: 4: 4= 4s 4: 4 j 4= 4t 4: 4: 4s 4: 4: 4: 4= 4s « >4 4t 4= 4t 4: 4: 4= 4= 4= 4: 4= 4e 4= 4< * 4t 4t 4: 4: 4s * 45 4= 4s 4s 4! 4t 4: 4« 4: 4s 4s 4= 4: :■ 


FOLLOWING USER'S GUIDES ARE PROVIDED BY DR. B.S.KIM AND ALL THE EQUATION 
NUMERS CORRESPOND EQUATIONS IN THE ACCOMPANYING TECHNICAL REPORT 


1. DATA TO BE SUPPLIED BY 4 DATA CARDS (IN THE ORDER OF READING) 


lANGLE: FIBER LAYUP ANGLE IN COMPOSITES 

C11,C12.. : ELASTIC MODULI OF COMPOSITE LAYERdN PSD 

NLAYER: NUMBER OF THE LAYERS IN THE GIVEN PLATE 
D, DELTA; THICKNESS OF COMPOSITE PLATEdN CM) 

RHO; MASS DENSITY OF COMPOSITE LAYERdN GR/CM«*3) 

NA; radius of THE IMPACT AS A MULTIPLE OF THE PLATE THICKNESS 
TAUO; impact TIME(IN SECOND) 

NX; NUMBER OF INEGRATION STEPS IN SPACE 
NT; NUMBER OF INTEGRATION STEPS IN TIME DOMAIN 
NXIMP; NUMBER OF SPACE STEPS IN IMPACT RADIUS 
NTIMP; NUMBER OF TIME STEPS IN IMPACT TIME 


2. WITH THE ABOVE DATA FOLLOWING PRIMARY DATA ARE CALCULATED 
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MAIN 


DATE = 77139 21/11/39 


C 

C 

C 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

r 

c 

c 

c 

c 

c 


3. 


4. 


5. 


B: A HALF OF THE LAYER THICKNESS 

K: wave number FOR FOURIER TRANSFORM 
S: LAPLACE TRANSFORM VARIABLE 

KO: LIMIT OF INTEGRATION FOR INVERSE FOURIER TRANSFORM FOR X 
OMEGAO: LIMIT 0"= INTEGRATION FOR INVERSE TRANSFORM IN TIME 
CO: LAPLACE TRANSFORM PARAMETER 


DISPLACEMENT AND STRESS FIELDS ARE CALCULATED IN TRANSFORMED SPACE 
AND BY INVERSIONS THESE BECOME DISPLACEMENTS AND STRESSES. 

THEY ARE GIVEN IN A MATIX FORM AS XX(I,J) REPRESENTING QUANTITY AT 
ITH TIME STEP AND JTH SPACE STEP IN X 

U(I,J): HRIZONTAL DISPLACEMENT 
V(I,J): VERTICAL DISPLACEMENT 
SIGMA(I,J): NORMAL STRESS 
TAU(I,J): SHEAR STRESS 
SIGMAKItJ): TANGENTIAL NORMAL STRESS 


FOLLOWINGS ARE WORKING MATRICES FOR THIS PROGRAM 


DATA(I,J), SUB(I,J): WORKING MATRICES FOR SUBROUTINE FOURT 
QO(I,J): INTEGRAL TRANSFORM OF IMPACT FUCTION GIVEN IN EQ(II-22) 

CAX(I,J), C8X(I,J): COS(ALPHA) AND COS(BETA) IN EQ(II-16) 
XlBXdfJ), .Y1AX( I T JJ , . . : XKBETA), YKALPHA),.. IN EQ(II-18) 


FOLLOWING SUBROUTINE ARE SUPPLIED IN THE PRESENT PROGLAM 

DPHASE: CALCULATES COS(BETA) AND CCS(ALPHA) IN EQ( 11-16) WITH GIV 
VALUES OF WAVE NUMBER K AND LAPLACE TRANSFORM VARIABLE 

dell: CALCULATES XKEETA), YKALPHA),.. IN EQ (11-19,20) 


ORIGINAL PAGE Ih 
DE POOR QUAUTY 
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DET; CALCULATE D,D1,.. IN EQ(II-25) 

FLIP: THIS PROGLAM GENERATES ONLY A HALF OF THE DATA (X>0) DUE 
TO SYMMETRY AND FLIP THEM TO FIND FULL DATA ALONG X AXIS 

MAP: CONTROLLS THE PRINTOUT FORMAT 

FACT: TAKES CARES OF THE COORDINATE SHIFT AND MULTIPLICATION FACTORS 
IN FOURIER-LAPLACE DOUBLE INTEGRAL TRANSFORMS 

FGURT: IS FAST FOURIER TRANSFORM ROUTINE SUPPLIED BY IBM 


C 

c 

IMPLICIT REAL*8(K) 

COMPLEX DATA I 32, 64) , SUB (32, 32) ,C AX ( 32 , 32 ) , CEX ( 32 , 32 ) , QO ( 32 ,32 ) 
COMPLEX S I GMA (32,32) ,SIGMA1( 32,32) ,TAU( 32,32) 

COMPLEX U( 32,32 ) ,V( 32,32 ) ,VX(32,32) 

COMPLEX Y1AX(32,32) , Y2 AX ( 32 , 32 ) , Y BAX ( 32 , 32 ) 

COMPLEX X13X( 32,32) , X2BX( 32 ,32 ) , X3BX ( 32 ,32) 

DIMENSION NN(2) ,MM(32) 

C 

REAL* 8 C11,C12,C22,C66,CHAT,E66, V66,PI,P2 
REAL* 8 DCOS,OSIN,DBLE,DSQRT,DFLOAT,DLOG 

REAL* 8 DELTA ,D,RH0,FN,F,KX(32) ,B,TAUO, A,DX,DT,UNITT,UNITX 
REAL*8 Q,OMEGAO,BL,CO 
C 

C0MPLEX*16 CDEX P, COLOG, CDSQRT,CDCOS, CDS IN 
C0MPLEX*16 BETA,ALPHA,CB,CA,Se,SA,C2NB,C2NA,S2NB,S2NA 
C0MPLEX*16 S,SI ,S2,D1,D2,D3,D4 
COMPLEX* 16 YIA, Y2A,Y3A,X18,X2B,X3B 
C 

COMMON YIA ,Y2A,Y3A,X1B,X2B,X3B,D1, D2,D3 ,D4,C11,C12,C22,C66,CHAT 
EQUIVALENCE (DELTA, D) 

EQUIVALENCE ( DATA ( 1 , 1 ), SUB ( 1 , 1 ) ) 

C 

SI=(0.D 00,1.D 00) 

PI=3. 1415926536D 00 
P2=PI*2.D 00 
INDEX=0 
INDEX=1 
C 

WRITE(6,4) 

4 FORMAT! • 1' ///////////////20X, •*«* WAVE PROPAGATION IN COMPOSITE P 

ILATE =!=**•) 
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MAIN 
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WRITE(6.5) 

5 =0RMAT(22X, 'GRAPHITE F I BER ( 55% )-EPCXY MATRIX COMPOSITE') 

C 

C jj: s;: j;: jj: j|: « sj! jji *«::{: :j: ^ ^ « 4! sji « 4: « 5): ♦=*=*« ’S' *=*=« =1= ^ #*****«=!'=!=* =i= * 

c 

c 

C INPUT data for elastic PROPERTIES OF COMPOSITE PLATE 
C ALL THE data ARE SUPPLIED IN PSI UNIT BUT NORMALIZED BY C66 
C WHICH IS CONSTANT REGARDLESS THE LAYUP ANGLE 
C 
C 

READ (5,101) IANGLE,C11,C12,C22,C66 

101 FORMAT! 110, 4D15. 7) 

CHAT=C 11-C12*«2/C22 

IF ( I ANGLE. EQ. 100) GO TO 200 
WRITE (6,102) IANGLE,C11,C12,C22,C66 

102 FORMAT!// 20X, 'LAYUP ANGLE = • , 1 3 , 3X , ' D EGR E E ' 

$ /20X, 'C( 1, 1)=' ,D12.5, ' PSI* tlOX, 'C(l,2)=' ,D12.5, ' PSI' 

$ /20X, 'C(2,2)=' ,D12.5, ' PS I • , lOX , ' C ( 6 , 6 ) = ' , D12 . 5 , ' PSI'//) 

GO TO 201 

200 CONTINUE 

WRIT£(6,210) Cll,C12,C22tC66 

210 F0RMAT(/20X, ' PLATE IS ISOTROPIC WITH POISSCN"S RATIO 1/4' 

S /20X, 'C(l,l) = ' »D12.5, • PSI'tlOX, 'C(l,2) = ' ,D12.5,' PSI' 

$ /20X, *C(2,2) = ' ,012.5, ' PS I • , lOX , • C (6 , 6 ) = • , D12 . 5 , • PSI'//) 

201 CONTINUE 
C 

E66=C66#6892.2D 00 
C11=C11/C66 
C12=C12/C66 
C22=C22/C66 
CHAT=CHAT/C66 
C66=l .D 00 
C 

C^f: WITH CORRECTION FACTOR 
C 

C66=P I««2/12.D 00 
C22=C22*C66 
C 
C 

C 

C 

C INPUT data for geometry OF COMPOSITE PLATE 

C ALL THE DATA ARE FIRST SUPPLIED IN CGS UNIT BUT CONVERTED INTO MKS UNIT 

C 

C 

READ(5,120) NLAYER, DELTA, RHO 
120 F0RMAT(I10,2020.10) 
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NL1=NLAYER+1 

FN=DFL0AT(NLAYER) 

B=DELTA/FN 

WRITE (6, 121 ) DELTA,RHO,NLAYER,e 

121 F0RMAT(20X, ' TOTAL THICKNESS OF COMPOSITE PLATE ; DELT A= • , F 10 . 5 , 

$ • CM'/ 

$ 20X, 'DENSITY OF COMPOSITE ; RHO= ' ,F10. 5, • GR/CM«>f'3'/ 

$ 20X, 'PLATE IS MADE OF ',I3,3Xt 'IDENTICAL LAYERS'/ 

$ 20X, 'LAYER THICKNESS ; 2B=',F10.5t' CM'//) 

C 

DELTA=DELTA/100.D 00 
RH0=RH0*1000.D 00 
B=B/200.D 00 
C 

c 

c 

c 

C INPUT DATA FOR IMPACT 

C 

C 

READ(5,60) NAtTAUO 

60 FORMAT!! 10tD20.10) 

READ(5,111) NX,NT,NXIMP,NTIMP 

111 F0RMAT(4I10) 

C 

WRITE(6,112) NX,NXIMP,NT,NTIMP 

112 F0RMATI20X,' TOTAL SPACE STEPS; NX = * , 1 3 , 5X , ' W ITH' , 1 3 1 2X , ' STEPS FOR 
S CONTACT RADIUS'/ 

$20X, 'TOTAL TIME STEPS ; NT= • , 13 » 5X , ' W I TH • , I 3 ,2X, ' ST E P S FORCONTAC 
$T TIME'//) 

C 

A=DFLOAT(NA)*DELTA 

DX=A/DFLOAT(NXIMP) 

DT=TAUO/DFLOAT(NTIMP) 

C 

WRITE(6,61) A,DX,TAUO,DT 

61 F0RMAT(20X, 'CONTACT RADIUS ; A='tD12.5,' M'/ 

$20X, 'SPACE STEP ; 0X = ',D12.5t' M'/ 

$20X, ' CONTACT TIME ; TAUO= • » D12 .5 , • SECOND'/ 

$20X,'TIME STEP ; DT=',012.5t' SECOND'//) 

C 

C 

C 

C 

C NORMALIZE ALL THE INPUT DATA 

C 

C 


V66=DSQRT{ E66/RH0) 
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UNITT=A/V66 
UNITX=DELTA 

F = D«D*CH0/ ( E66*UNITT=f=«2 )/2.D 00/ FN 
C 

NX2=NX/2 
NN( 1 )=NT 
NN( 2) =NX 
DX=DX/UNITX 
DT=DT/UNITT 
K0=PI/DX 
OMEGAO=PI/DT 
BL=A/D 

C0=0L0G(l-0 06=f=2.D 00*DT)/(3.D 00*DT=i=DFL0AT( NT) ) 






CALCULATE THE IMPACT INPUT FUNCTION Q0(I,J) IN EQ(II-22) 

CALCULATE C0S(6£TA) AND COS( ALPHA) IN EQ( 11-16) BY SUBROUTINE DPHASE 
CALCULATE Xl(BETA), YKALPHA) BY SUBROUTINE DELL 


DO 30 J=i,NX2 

K=2.D 00*K0*(DFL0AT( J )-.5)/0FL0AT(NX)-K0 

K2=K**2 

KX( J)=K 

Q = PI«*2/DSQRT(P2)*DSIN(K*BL)/K/( ( K*BL ) **2-P I =«=*2 ) 

C 

DO 30 1=1, NT 

S = C0 + SI«0MEGA0*(1.0 00- ( DFLOAT ( I ) - . 5D 00)«2.D OO/DFLO AT I NT ) ) 
S2=S**2«F 

QO(I , J ) = Q/2./S*( 1 .D 00-CDEXP(-S-TAU0/UNITT) ) # ( PB’^UNI TT ) 

$ /( ( S=»=TAU0)**2+(P2*UNITT)=!=*2) 

C 

CALL DPHASE ( K , S2 , C B , C A , NLAYER ) 

CBX( I , J)=C8 
CAX(I,J)=CA 
C 

CALL DELL(K,S2,CB,CA,SI,NLAYEP ) 

Y1AX( I , J )=Y1A 
Y2AXI I ,J)=Y2A 
Y3AX( I ,J)=Y3A 
X1BX( I , J)=X1B 
X2BX( I ,J)=X2B 
X3BX( I , J)=X3B 
C 
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30 CONTINUE 
C 
C 
C 

c 

c 

C REPRODUCE THE IMPACT FUNCTION TO CHECK INPUT 
C 

DO 300 1=1 ,NT 
DO 300 J=1,NX2 

300 SUB( I , J) =Q0( I » J ) 

CALL FLIP{DATA,NX,NX2,NTt+l) 

CALL F0URT(DATA,NN,2,-ltl,0) 

CALL FACT(DATA,NX,NT,CO,OMEGAO,KO,PI,SI ) 

WRITE! 6,301) 

301 FQRMATC l'///////////////20X, REPRODUCTION OF IMPACT FUNCTIO 

IN «*«*) 

CALL MAP(DATA,NX,NT,NX2 ,MM,INDEX) 

C 

C 

c 

c 

c 

C THIS IS THE MAIN PART OF THE PROGRAM. 

c 

C CALCULATE D( BET A ) , OBAR ( AL PHA ) , . . IN EQ(II-19,20) BY SUBROUTINE DET 
C NEXT CALCULATE 1,V,.. IN EQ(II-18) IN TRANSFORMED SPACE 

C AND FLIP TO FIND FULL DATA AND INVERT THEM BY MEANS OF FOURT. 

C REPEAT THIS PROCESS FROM N=0 TO NLAYER 
C 

DO 11 N=1,NL1 
NY=N-1 
NYY=NY-1 
C 
C 

C 

C 

C GENERATION OF DATA FOR DISPLACEMENTS AND STRESS IN TRANSFORMED SPACE 
C 

DO 100 J=1,NX2 
DO 100 1=1, NT 
CB=CBX( I ,J) 

CA=CAX( I ,J) 

X1B=X1BX(I ,J> 

X2B = X2BXI I ,J ) 

X3B=X3BX( I ,J ) 
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Y1A=Y1AX( I ,J) 

Y2A=Y2AX( I , J) 

Y3A = Y3AXn ,J) 

C 

SB=CDSQPT( l.D 00-CB**2) 

SA=CDSGRT(1.0 00-CA**2) 

BETA=C6+SI*SB 
ALPHA=CA+SI*SA 
BETA=CDLOG{BETA)/SI 
ALPHA =CDLOG( ALPHA) /SI 
C 

CALL 0ET(ALPHA,8ETA,SI,FN) 

C2N3=CDCCS (2.D 00*8 ETA«DF LOAT ( NY ) ) 

S2N3 = CDSIN(2.D 00*B ET A=('DFLOAT ( NY ) ) 

C2NA = CDCOS(2.D 00«ALPHA=f=DFLCAT ( NY ) ) 

S2NA = CDSIN(2.D 00* ALPHA=(=DFLOAT (NY ) ) 

C 

U( I, J ) = (X1B*{01*C2NB + SI*D2*S2NB)+Y1A*(04*C2NA+SI*D3*S2NA) )*Q0( I ,J) 
V( I , J) =(X2B*(D2*C2NB+SI*D1*S2NB)+Y2A*(D3*C2NA+SI*D4*S2NA) )*Q0( I,J) 
TAU( I ,J)=( X3B*{D2*C2NB + SI*Dl*S2NB) + ( D2*C2NA+SI*D4*S2NA) )*Q0( I, J) 
SIGMA( ItJ)=( (D1*C2NB+S I*D2*S2NB)+Y3A*(D4*C2NA+SI*D3*S2NA) )*Q0( I,J) 
100 CONTINUE 
C 
C 

C 

C INVERSION AND PRINTOUT OF HORIZONTAL DISPLACEMENT UN(ItJ) 

C 

DO 10 1=1, NT 
DO 10 J=1,NX2 
10 SUB( I , J)=U( I ,J ) 

CALL FLIP(0ATA,NX,NX2,NT,-1) 

CALL FOURT(DATA,NN, 2,-1, 1,0) 

CALL FACT! DAT A , NX , NT , C 0 , OMEGAO , KO , P I , S I ) 

WRITE(6,981) NY 

981 F ORM AT {•!•// ////////// ///20X,*U* ,13) 

CALL MAP (DATA, NX, NT ,NX2,MM, INDEX) 

C 

C 

C 

C INVERSION AND PRINTOUT OF VERTICAL DISPLACEMENT VN ( I , J ) 

C 

DO 20 1=1, NT 
DO 20 J=1,NX2 
SUB( I , J)=V( I , J) 

CALL FLIP(DATA,NX ,NX2,NT,+1) 

CALL FOURT(DATA,NN, 2,-1, 1,0) 


20 
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CALL FACTI DATA,NX ,NT,C0,0MEGA0,K0, PI ,SI ) 

WRITE{6t982) NY 

982 FORMAT! • 1* /////// //////// 20X, • V ,13) 

CALL MAP(DATA,NX,NT,NX2, MM, INDEX) 

C 

C 

C 

C INVERSION AND PRINTOUT OF SHEAR STRESS TAU(I,J) 

C 

DO 35 1=1, NT 
DO 35 J=1,NX2 
35 SUB{ I , J)=TAU( I , J) 

CALL FLI P(DATA,NX ,NX2,NT,-1) 

CALL FOURT(DATA,NN, 2,-1, 1,0) 

CALL FACT(DATA,NX,NT,CO,OMEGAO,KO,PI,SI ) 

WRITE(6,983) NY 

983 FORMAT!' 1 '///////// //////20X ,* TAU • ,13) 

CALL MAP !DATA,NX,NT ,NX2,MM, INDEX) 

C 

C 

c 

C INVERSION AND PRINTOUT OF NORMAL STRESS SIGMA!I,J) 

C 

DO 40 1=1, NT 
DO 40 J=1,NX2 
40 SUB! I ,J)=SIGMA! I,J) 

CALL FLIP!DATA,NX,NX2,NT,+1) 

CALL FCURT!DATA,NN,2,-1,1,0) 

CALL FACT!DATA,NX,NT,CO,OMEGAOrKOrPI ,SI ) 

WRITE!6,984) MY 

984 FORMAT!' 1 '//////////// ///20X ,' S I GMA ' ,13) 

CALL MAP! DAT A, NX, NT ,NX2 , MM , I NDEX ) 

C 

C 

C 

C INVERSION AND PRINTOUT OF TANGENTIAL NORMAL STRESS SIGMA1!I,J) 

C 

IP INY.EQ.O) GO TO 160 
DO 50 1=1, NT 
DO 50 J=1,NX2 

50 SUB! I , J)=SIGMA1! I , J ) + FN*C12*V! I, J) 

CALL FLIP!0ATA,NX,NX2,NT,+1) 

CALL FOURT!OATA,NN,2,-1,1,0) 

CALL FACT!DATA,NX,NT,CO,OMEGAO,KO,PI,SI ) 

WRITE!6,985) NYY 
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985 FORMAT! • 1' //////////// ///20X, • SI GMAl ' , 13) 

CALL MAP(DATA,NX,NT,NX2,MM,INDEX) 

C 

IF ( NY.EQ.NLAYER) GO TO 70 
DO 51 1=1, NT 
DO 51 J=1,NX2 

51 SIGMAl ( I ,J )=-SI«KX(J) «C11-U( I , J)-FN=!'C12=(=VX( I , J) 

GO TO 80 
C 

160 CONTINUE 

DO 161 1=1, NT 
DO 161 J=1,NT 

161 SIGMAK I ,J)=-SI*KXI J)=!'C11*U( I , J )-FN=i=C12«V { I , J ) 

GO TO 80 

C 

70 CONTINUE 

DO 71 1=1, NT 
DO 71 J=1,NX2 

71 SUB( I , J) =-SI*KX( J )*C11«U( I , J) + FN*C12=«=( V( I , J)-VX( I , J) ) 

CALL FLIP(DATA,NX,NX2,NT,+1) 

CALL F0URT(0ATA,NN,2,-1,1,0) 

CALL FACT! DATA, NX,NT,C0,0MEGA0,K0,PI ,SI ) 

WRITE!6,985) NY 

CALL MAP(DATA,NX,NT,NX2,MM, INDEX) 

GO TO 90 
C 

80 CONTINUE 

DO 81 1=1, NT 
DO 81 J=1,NX2 

81 VXI I , J)=V! I, J) 

C 

90 CONTINUE 

C 

C 

11 CONTINUE 

C 

^ ^ ^ jJ: :(e 4= 5|: is * ^ : 

C 

C 

c 

c 

STOP 

END 
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#:!£ jjs !j! :j! 4! * =i: ^ ###*#*« sf:**** sjs =i: **=!=** =r ^ 

THIS SUBROUTINE CALCULATES THE PHASE SHIFT BETA AND ALPHA FROM 

EQdI-15,16) OF THE PRESENT REPORT WITH GIVEN VALUES OF 
WAVE NUMBER K AND LAPLACE TRANSFORM VARIABLE S 

CA. = COS( ALPHA) 

CB=COS(BETA) 

jj: 4: :(s sj: sj: 3): Jj: Jjt « ^ 5|t 5f! J|! # X: >f= >r *1= sf: ^ ***«# ^: **=!= 5(: * 


SUBROUTINE DPHASE(K,S , CB ♦ C A,NLAY EP ) 

IMPLICIT C0MPLEX*16 (A, XtY) 

COMPLEX* 16 ROOTP ,ROOTM,StCDSQRT,DCMPLX ,DCB,DCA 
C0MPLEX*16 DlfD2f03,D4 
C0MPLEX*16 CB,CA 

REAL* 8 C11,C12,C22,C66,CHAT,N,K2,DFL0AT,DBLE,K 
C 

COMMON YlA,Y2A,Y3A,XlB,X2BtX3B,Dl,D2,D3,D4,Cll»C12,C22,C66,CHAT 
ROOTPI AA,ABtAC)={-AB+CDSQRT{AE**2-4.D 00*AA*AC) ) / ( 2.D 00*AA) 
ROOTMI AA, AB, AC) =(-AB-CDSQRT(AB**2-4.D 00*AA*AC ) ) / ( 2. D 00*AA) 

C 

N=DFL0AT(NLAYER)*2.D 00 
K2=K**2 
C 
C 

A1=(K2*C11/N+S ) *(K2*C66/N+S) 

A2=K2*CHAT/(3.D 00*N)+N*C66+S/3. D 00-C 12*C66*K2/ ( 3 . D 00*N*C22) 
A2=A2*(-C12*K2/ (3.D OC*N) +N*C22+S/ 3. D 00) 

A3=N*C22+S/3.D 00-K2*C12* (C66*K2/N+S )/ ( 9. D 00*N**2*C2 2 ) +C66*K2/ 
$ (3.D 00*NJ 

A3=A3*(C11*K2/N+S)-K2*(C12+C66)**2 

A3=A3+(C66*K2/N+S)*{CHAT*K2/(3.D 00*N) +N*C66+S/3 . +C 12 **2*K2/ 

$ (3.D 00*C22)) 

C 

AA=A1+A2-A3 
AB=A3-2.D 00*A2 
AC = A2 

DCB=ROOTPIAA,AB»AC) 

DCA=ROOTM( AA,AB,AC) 

C 

CB=CDSQRT(DCB) 

CA=CDSQRTI DCA) 

RETURN 
END 


ORIGINAi; PAGE IS 
OF POOR QUALITYI 
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G1 RELE/^SE 2.0 


MAIN 


DATE = 77139 21/11/39 


« jj: J|! Jjssit « 4! ii: « #3!: si! jft «*«**««# si: s;: ** *5!^ «««**««** ^ * 

THIS SUBROUTINE CALCULATES DELTA ( BET A) , DELTABAR ( ALPHA ), DELTAl ( BETA ),. . 
IN EQ(II-19,20) AND XI ( BETA ) , Y 1 ( ALP HA ) IN EQ 1 1 I -1 8 , 19 , 2 0) 

X1B=X1(BETA) » Y1A=Y1( ALPHA) 


X= 3 «: 3«: =!3 3f: =!3 3f! * =S= =!= 3|: « 3(3 «««#«#« =1= * 31: Jis** >!5 « 


SUBROUTINE DELL ( K , S , C8 , C A , SI » NLAY ER ) 

IMPLICIT REALT-8(K) , COMPLEX^ 16{ S , X , D , Y ) 

C0MPLEX3(316 SB,CB,CDSQRT, SA,CA 
REAL*8 N,DFLOAT 
REALMS C11,C12»C22,C66,CHAT 
C 

COMMON YlA,Y2A,Y3A,XlE,X2B,X36,Dl,D2»D3,D4,ClltC12»C22,C66,CHAT 
C 

K2=K«*2 

N=DFLOAT(NLAYER) 

C 

C 

Sll = S+Cll>i=K2/2.0 00/N 
S66=S+C66*K2/2.D 00/N 
SA33=CDSGRT(1.D 00-CA*«2) 

SB=CDSQRT( l.D 00-C5«*2) 

C 

C 

DELTAB=CB*’!=3*S1 l*S66+SB<=«2=f=CB*{-C66«lC66 + C12 J*K2 
$ +S66*(S/3.D 00+CHAT*K2/6.D 00/N+2.D 00=!'N=!'C66 ) ) 

DELie=SI*K=«=SB**2’!=CB=!'(-C66-C12+C12«S66/6.0 00/N/C22 ) 

DEL2B = SI’!=S3**3*(C12*C66*K2/6.D 00/ N/C22-S/3 . D 00 
$ -CHAT«K2/6.D 00/N-2.D 00*N*C66 )-S I >!=CB*=!=2*S B*S 1 1 
DEL3B=CB«*2«SB«(C12*K*S11*S66/6.D 00 /N/C2 2-C 1 1 ) 

$ +SB**3*^C12«K*{ S/3.D 00+CHAT*K2/6.D 00/N+2.D 00*N*C66) 

$ -C12**2*C66*K*K2/6.D 00/N/C22) 

C 

0ELTAA=3:SI=!3CA«*2*SA*( ( C12 + C66)*C66*K2-S66=!'( S/2.D 00+CHAT=!=K2 
$ /6.D 00/N + 2.0 03*N*C66+C12’«=«2*K2/6.D 00/N/C22)) 

$ +Si*SA=i==(=3*(S/3.D 00+2. D 00*N«C22 ) * ( C66 *C 1 2-K2/6 . D 00/N/C22 

$ -S/3.D 00-CHAT’»=K2/6.D 00/N-2.D OO’i^N^Cdb) 

DEL1A=SA*>!»2«CA=(=(-K2*C12*S66/6.D 00/6. C 00/N’i=*2/C22+C66*K2/6. D 00/N 
$ +S/3.D 00+2. D 00«N*C22)+CA«*3*S66 

DEL2A = CA*4=2*SA=J'K*(C12+C66) 
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G1 RELEASE 


C 


2.0 DELL DATE = 77139 

+ S A*=«'3’!'(K/6.D 00/N=«' ( S/3.D 00+CHAT’f=K2 /6 . D 00/N+2.D 
-C66*Cl2*K=f'’!'3/36.D 00/C22/N>i'«2 ) 
DEL3A=-SI*CA’l=*3=S'C12*K’«=S66 + SI*SA--2*CA*( C66-K*( S/3.D 0 
*C22+C66>i‘K2/6.D 00/N)-K/6.D 00/N«S66«( S/3 .D 00+CHA 
/N+2.D 00*N*C66) ) 


X1B=-DEL1B/DELTAB 

X2B=-DEL2B/DELTAB 

X3B=-DEL3B/DELTA6 

Y1A=-DEL1A/DELTAA 

Y2A=-DEL2A/DELTAA 

Y3A=-DEL3A/DELTAA 

RETURN 

END 


21/11/39 

00«N*C66) 

0+2. D 00=«'N 
T=!=K2/6.D 00 
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RELEASE 2.0 


MAIN 


DATE 


77139 


21/11/39 


3^!-: 5^5 «sjt« * 5 ;: jj:« 3!: 3js3e:3j:3!:3i!3{:3i:3f[3|!3i:3!:3i£3!:3!c3!::S:s;:3i:3js3!::S::jt3f:#!«:3!!3i!3j::(s3!3 

THIS SUBROUTINE CALCUALTES D»D1,.. IN EQ( 11-23) OF THE PRESENT REPORT 

0 « 3!: « 3|C ^ 3ic 3h !«: 3!t 3(! 3j! 3[: * 3js !(: 3j! 3j: 3>£ «« 3!t :t: # 3(: 3(3 3|t ««« 8(: ««« 45 «*#«««« =i: «*=«:« Jf: ^ sj: 3>3 * 3js :{: 5!: 3it 3!: « 3is 3js 

C 

C 

SUBROUTINE DET { ALPH A , B ET A , S I , FN ) 

IMPLICIT C0MPLEX*16(D»XtY] 

C0MPLEX*16 ALPHA, BETA 

COMPLEX* 16 C2N6,C2NA,S2NA,S2NB,CDSQRT, SI 
REAL*8 FN 

REAL*8 C11,C12,C22,C66,CHAT 
C 

COMMON YIA, Y2A,Y3A,X16,X2B,X3B,D1,D2,D3,D4,C11,C12,C22,C66 tCHAT 
C 

C2NA=CDC0S(2.D 00*ALPHA*FN) 

S2NA5CDSIN(2.D 00*ALPHA*FNJ 
C2NB=CDC0S(2.0 00*BETA*FN) 

S2NB=CDSIN(2.D 00*BETA*FN) 

C 

X=Y3A*X3B*(1.D 00-C2NA*C2NB ) 

Y=X3B*Y3A*S2NB*C2NA-S2NA*C2NB 

C 

D=-2.D 00*X+(1.D 00+X3B**2*Y3A**2)*S2NA*S2NB 
Dl=- ( X-S2NA*S2NB) 

D2=-SI*Y 

D3=SI*Y*X3B 

D4=X3B*( X3B*Y3A*S2NB*S2NA+C2NA*C2NB-1.D 00) 

C 

D1=DI/D 

02=02/0 

03=03/0 

04=04/0 

C 

RETURN 

END 
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G1 PELEASE 2.0 MAIN DATE = 77139 21/11/39 

C 

C 


ALL THE DATA IN THE MAIN PROGRAM ARE GENERATED FOR ONLY HALF OF THE 
PLATE WHEN X>0 . DUE TO SYMMETRY OF THE PROBLEM WE CAN GENERATE 
THE FULL DATA BY FLIPPING THE HALF OF THE DATA. 


:(= 3fr « 3(! :):* :j:5S!3i::S: a): « 3J: ««#*#««««« «« 5!: Jj! « ^ 5i! 5|s « *:(£ :^s 


SUBROUTINE F L I P ( DATA , NX , NX2 t NT , I NDEX ) 
COMPLEX DATA(NT,NX) 

DO 10 J=1,NX2 
JJ=NX+1-J 
DO 10 1=1, NT 

DATA! I,JJ)=FLOAT{ I NDEX ) *DATA ( I , J ) 

10 CONTINUE 
RETURN 
END 
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G1 RELEASE 2.0 MAIN DATE = 77139 21/11/39 

C 

C 

jj: sic :(! ^5(1 4= sjc:;: sjcsjc s^si: * sic sfcsjcsfc*:}: *si: sjcsjc sjisfcsjc sj: 

c 

c 

C THIS SUBROUTINE TAKES CARE OF THE COORIDINATE SHIFT IN 
C LAPLACE AND FOURIER TRANSFORM IN THE PROCESS OF APPLYING 

C FAST FOURIER TRANSFORM ALGORITHM 

C 
C 

C sjc s}: sic s(c !{: s|c * s{c * s{c « * 3(C s<c « 5|C sjc s(c jjc sic « s(s S{C sit # sic s(c Stc sjc sjc « « 3{C « ijc s(c * # S|C sic « sic « s(c « !(C B(C « SfC sic sic s|c sj: !(C S}! sic sic sic « :«c « sic sic s!c 3|C sic sic sic :jc S{C sic 

c 

c 

SUBROUTINE F AC T ( OAT A , NX , NT , CO » WO , KO , P I » S I ) 

COMPLEX DATA(NTtNX) 

C0MPLEX*16 CDEXP,SI 
REAL*8 DEXP,DSQRT,DFLOAT 
REAL*8 COtWO,PI ,KO tFT,FX 
FX=DFLOAT( NX) 

FT=DFLOAT(NT) 

NX2=NX/2 

C 

DO 10 L=ltNX2 
DO 10 M=1»NT 

0ATA(M,L)=DATA(M,LJ*4.D 00*KO*WO/ { DSQRT ( 2. D 00*PI )«*3*FT«FX) 

$ *DEXP(C0*PI#DFL0AT(M-1)/W0)«CDEXP(SI«PI«{ l.D 00-1. D 00/FX) 

$ *DFLOATlL-in*CDEXP(SI«PI«{l.D 00-1. D 00/FT)*DFLOAT( M- 1 ) ) 

10 CONTINUE 
RETURN 
END 




OY 
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Gl RELEASE 2.0 main DATE = 77L39 21/11/39 

C 

C 

C 

C THIS SUBROUTINE CONTROLS THE FORMAT OF THE PRINTOUT OF THE FINAL RESULTS 
C 

C IF INDEX=0: ALL THE NUMERICAL VALUES ARE PRINTED 

C =1; THE MAXIMUM AND NORMALIZED VALUES ARE PRINTED 

C 
C 

s}: !{t * >!t :Je « :ic :{s Jjt :it :(s jjt j<: jj! ^ jjs j|t ^ :{t a}: jjt :js ){: «!(: <s Jj: #s{:* 5«! :ji « :J: # :j£ 4s 

C 

C 

SUBROUTINE M AP ( D AT A , NX , NT ,NX2» MM , I NDEX ) 

COMPLEX DATA(NT,NX) ,S 
DIMENSION MM(NX2) 

IF ( INDEX. EQ.l) GO TO 200 
DO 44 IQ=1,NT 

44 WRITE(6tl5) I Q , ( D AT A ( I Q , I) , I =1 ,NX2 ) 

15 FORMAT! 15, 2E 14.5 t2X,2E 14.5, 2X,2E 14.5, 2X,2E 14. 5/ 

$ 3(5X,2E14. 5,2X,2E14.5,2X,2E14.5,2X,2E14.5/)/ 

$ 4(5X,2E14. 5,2X,2E14.5,2X,2E14.5,2X,2E14.5/)/) 

200 CONTINUE 

find THE MAXIMUM VALUE 

RS= l.E-3 

NT5=NT-5 

NX5=NX2-5 

DO 114 1=1, NT5 

DO 114 J= 1,NX5 

S= DATA! I, J) 

TP= PEAL(S)/RS 

IF (ABSITPJ.LT.l.) GO TO 114 
RS= REAL(S) 

114 CONTINUE 

WRITE (6,516) RS 

516 FORMAT! 20X, •=!=** MAXIMUM VALUE =',E12.5,' «4=*»//) 

DO 119 1=1, NT5 
DO 113 J=1,NX5 
S= DATA! I, J) 

113 MM!J)= REAL! S) /PS*100 

WRITE!6,515) ( MM ! K I M) , K IM=1 , 27 ) 

119 CONTINUE 

515 FORMAT! lOX, 2713) 

RETURN 

END 
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1 RELEASE 2.0 FOUPT date = 771^9 

subroutine four T( data, NN t no IM, ISIGN, I form, WORK) 
DIMENSION DATA! 1) ,NN( 1) , IFACT( 32) ,WORK( 1) 

TWOPI =6.283185307 
IF(NDIM-1)920, 1, 1 
1 NT0T=2 

DO 2 IDIM=1,NDIM 
IF(NN( IDIM) ) 920,923,2 
NTOT=NTOT*NN( IDIM) 

MAIN LOOP FOR EACH DIMENSION 

NP1=2 

DO 910 IDIM=1,NCIM 
N=NN( IDIM) 

NP2=MP1*N 
IF(N-1 )920,900, 5 


21/11/39 


FACTOR N 


M=N 

NTWQ=NP1 
IF = 1 
IDIV=2 
/O IQUOT=M/IDIV 

IREM=M-IOIV*IQUOT 
IF( IQUOT-IDIV)50, 11,11 
II IF( IREM)20,12,20 
/X NTWO=NTWO+NTWO 
M=IQUOT 
GO TO 10 
2(7IDIV = 3 
JO IQUOT=M/IDIV 

IREM = M-IOI V’i^IQUOT 
IF( IQUCT-IDIV)60,31,31 
31 I^I IREM)4C,32,40 
31l FACT { IF) = IDIV 
IF=IF+1 
M=IQUOT 
GO TO 30 
^IDIV=IDlV+2 
_ GO TO 30 
tsOlFI IREM)60,51,60 
51 NTWO = NTWO+NTWO 
GO TO 70 
^IFACTI IFJ=M 

A 

C separate four CASES — 

Cj 1. COMPLEX TRANSFORM OR REAL TRANSFORM FOR THE 4TH, 5TH,ETC. 


FFTTOOO 

FFTT077C 

FFTT078.'- 

FFTT079v 

FFTT0801 

FFTT081 

FFTT082C 

FFTT083f, 

FFTT084 . 

FFTT085C 

FFTT086C 

FFTT087'. 

FFTT083C 

FFTT039": 

FFTT090C 

FFTT091C 

FFTT092 : 

FFTT093C 

FFTT094C 

FFTT095 . 

FFTT096C 

FFTT097C 

FFTT093C 

FFTT099C 

FFTTIOOC 

FFTTIOIC 

FFTT102C 

FFTT103I 

FFTT104C 

FFTTlOSi 

FFTT106v 

FFTT107C 

FFTT103C 

FFTT109C 

FFTTllO . 

FFTTlllC 

FFTT112C 

FFTT113: 

FFTT114C 

FFTT115C 

FFTTl 16C 

FFTT117C 

FFTT118 : 

FFTT119C 

FFTT120C 

FFTTi 210 

FFTT122C 

FFTT123C 
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1 RELEASE 2. 0 


FOURT 


DATE = 77139 21/11/39 


C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

70 


71 

72 

73 

74 


80 

90 

95 

C 

C 

C 

C 

100 

110 


120 


DIMENSIONS. FFTT124C 

2. REAL TRANSFORM FOR THE 2ND OR 3RD DIMENSION. METHOD — FFTT125C 

TRANSFORM HALF THE DATA, SUPPLYING THE OTHER HALF BY CON- FFTT126C 

JUGATE SYMMETRY. FFTT127C 

3. REAL TRANSFORM FOR THE 1ST DIMENSION, N ODD. METHOD— ' FFTT128C 

TRANSFORM HALF THE DATA AT EACH STAGE, SUPPLYING THE OTHER FFTT129C 
HALF BY CONJUGATE SYMMETRY. FFTT130C 

4. REAL TRANSFORM FOR THE 1ST DIMENSION, N EVEN. METHOD — FFTT131C 

TRANSFORM A COMPLEX ARRAY OF LENGTH N/2 WHOSE REAL PARTS FFTT132C 
ARE THE EVEN NUMBERED REAL VALUES AND WHOSE IMAGINARY PARTS FFTT133G 
ARE THE ODD NUMBERED REAL VALUES. SEPARATE AND SUPPLY FFTT134C 

THE SECOND HALF BY CONJUGATE SYMMETRY. FFTT135C 

FFTT136C 


N0N2=NP1*( NP2/NTW0) 

ICASE=1 

IF( IDlM-4) 71,90,90 
IF( IF0RMJ72, 72,90 
ICASE=2 

IF(IDIM-1)73,73,90 

ICASE=3 

IF(NTWC-NP1)90,90,74 

ICASE«4 

NTW0=NTW0/2 

N=N/2 

NP2=NP2/2 

NT0T=NT0T/2 

1 = 3 

DO 80 J=2,NT0T 
DATA! J)=DATA(I ) 

1 = 1+2 
I1RNG=NP1 

IF( ICASE-2)100,95,100 
I1RNG = NP0>!'( l+NPREV/2) 

SHUFFLE ON THE FACTORS OF TWO IN N. AS THE SHUFFLING 

CAN BE DONE BY SIMPLE INTERCHANGE, NO WORKING ARRAY IS NEEDED 

IF{ NT WO-NPl) 600,600,110 
NP2HF=NP2/2 
J = 1 

DO 150 I2=1,NP2,N0N2 
IF( J-I2) 120, 130,130 
IlMAX=I2+N0N2-2 
DO 125 11=12, I 1MAX,2 
DO 125 13=11 ,NT0T,NP2 
J3=J+I3-I2 
TEMPR=DATAC I3J 
TEMPI=DATA( 13+1 ) 


FFTT137C 
FFTT138C 
FFTT139C 
FFTT140C 
FFTT141C 
FFTT142C 
FFTT143C 
FFTT144: 
FFTT145C 
FFTT146C 
FFTT147C 
FFTT148C 
FFTT149Q 
FFTT150C 
FFTT151C 
FFTT1520 
FFTT153C 
FFTT154C 
FFTT1550 
FFTT156C 
FFTT157C 
FFTT158C 
FFTT159C 
FFTT160C 
FFTT161C 
FFTT162C 
FFTT1630 
FFTT164C 
FFTT165C 
FFTT166C 
FFTT167C 
FFTT168C 
FFTTi d90 
FFTT170C 
FFTT171C 
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, date 

= 77129 

21/11/39 


0ATA( I3)=DATA{ J3) 



«=FTT172 


DATA( 13+1) =DATA(J3+1) 



FFTT173 . 


DATA( J3) =TEMPR 



FFTT174: 

125 

DATA( J3+1) =TEMPI 



FFTT175 - 

130 

M=NP2HF 



FFTT176;; 

140 

IF(J-M)150t150,145 



FFTT1771 

145 

J = J-M 



FFTT173C 


M=M/2 



FFTT179C 


IF(M-N0N2)150,140,140 



FFTT180C 

150 

J=J + M 



FFTT181C 

C 




FFTT182C 

C 

MAIN LOOP FOR FACTORS OF TWO. 

PERFORM 

FOURIER 

TRANSFORMS OF FFTT183r 

C 

LENGTH FOUR, WITH ONE OF LENGTH 

TWO IF 

NEEDED. 

THE TWIDDLE FA CTOR FFTT 1 84- 

C 

W=EXP( ISIGN«2’!'PI«SQRT(-L) «M/ (4* 

MMAX) ) . 

CHECK 

FOR W=ISIGN«SQRT(-l)FFTT185v 

C 

AND REPEAT FOR W= I S I GN*SQRT ( -1 ) *CON J UG AT E ( W ) . 

FFTT;, 86 .. 

C 




FFTT187C 


N0N2T=N0N2+N0N2 



FFTT188C 


IPAR=NTW0/NP1 



FFTT189 ; 

310 

IF( IPAR-2)350,330,320 



FPTT190t 

320 

IPAR=IPAR/4 



FFTT191C 


GO TO 210 



FFTT 192: 

330 

DO 340 1 1 = 1, I1RNG,2 



FFTT193C 


DO 340 J3=I1,N0N2,NP1 



FFTT194C 


DO 340 K1=J3,NTOT,NON2T 



FFTT195C 


K2=K1+N0N2 



FFTT196C 


TEMPR=DATA(K2) 



FFTT 197: 


TEMPI=CATA(K2+1 ) 



FFTT198C 


DATAIK2) =DATA(K1)-TEMPR 



FFTT199C 


DATA( K2+1)=0ATA(K1+1)-TEMPI 



FFTT200C 


DATA( Kl) =DATA(K1)+TEMPR 



FFTT201C 

340 

DATA( K1+1)=DATA{K1+1)+TEMPI 



FFTT202C 

350 

MMAX=NCN2 



FFTT203-: 

360 

I F(MMAX-NP2HF) 370 ,600,600 



FCTT204' 

370 

LMAX=MAXO(NON2T,MMAX/2) 



FFTT205' 


I F(MMAX-NON 2)405,405,380 



FFTT206- 

380 

THETA=-TWOPI*FLGAT( N0N2)/ FLOAT (4*MMAX) 


FFTT207C 


IFI ISIGN)400,390,390 



FFTT 2 08: 

390 

theta=-theta 



FFTT209C 

40C 

WR=COS(THETA) 



FFTT210C 


WI=SIN(THETA) 



FFTT211: 


WSTPR=-2.*WI*WI 



FFTT212C 


WSTPI =2.*WR’!'WI 



FFTT213C 

405 

DO 570 L=N0N2,LMAX,N0N2T 



FFTT214C 


M=L 



FFTT215C 


IF(MMAX-N0N2) 420,420,410 



FFTT216C 

410 

W2R=WR*WR-WI*WI 



FFTT217C 


W2I = 2 .>!'WR«WI 



FFTT218C 


W3R = W2P«WP.-W2I*WI 



FFTT219C 
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W3I=W2R*WI+W2I*WR 
420 DO 530 I1=1,I1RNG,2 

DO 530 J3=IltN0N2,NPl 
K>IIN= J3+IPAR«M 
I F(MMAX-NDN2) 43 0,430,440 
430 KMIN=J3 
440 KDI>== IPAR^MMAX 
450 KSTEP=4«KDIF 

DO 520 K1=KMIN,NT0T,KSTEP 

K2=K1+KDIF 

K3=K2+KDIF 

K4=K3+KDIF 

IF{MMAX-N0N2) 460, 460,480 
460 U1R=DATA(K1)+DATA(K2) 

U1I=DATA(K1+1)+DATA(K2+1J 
U2R=DATA(K3J+DATA(K4) 
U2I=DATA(K3+1)+DATA(K4+1) 
U3R=0ATA(K1)-DATA(K2) 
U3I=DATA(K1+1)-DATA( K2+1) 

IF( ISIGN)470,475,475 
470 U4R=DATA(K3+1)-0ATA(K4+1) 

U4I=DATA(K4)-0ATA(K3) 

GO TO 510 

475 U4R=DATA(K4+1)-DATA(K3+1) 

U4I=DATA(K3)-DATA(K4) 

GO TO 510 

480 T2R=W2P«DATA(K2)-W2I*DATA(K2+1) 
T2I = W2R*0ATA(K2+l)-»-W2I*DATA(K2) 
T3R = WR>!=DATA{ K3 )-W I *DAT A ( K3+ 1 ) 
T3I = WR’«'DATA( K3+1 )+WI*DATA(K3i 
T4R=W3R«DATA(K4)-W3I«DATA{K4+11 
T4I=W3R«DATA(K4+1 )+W3I*DATA(K4) 
U1R = DATA( KD+T2R 
U1I=DATA(K1+1)+T2I 
U2R=T3R+T4R 
U2I=T3I+T4I 
U3R=DATA(K1)-T2R 
U3I=DATA(K1+1)-T2I 
IF( ISIGN)490,500,500 
490 U4R=T3I-T4I 
U4I=T4P-T3R 
GO TO 510 
500 U4R=T4I-T3I 
U4I=T3R-T4R 

510 DATA( K1)=U1R+U2R 

DATA(K1+1)=U1I+U2I 
DATA( K2)=U3R+U4R 
DATA( K2+1) =U3I+U4I 


DATE = 77139 21/11/39 

FFTT223 

FFTT221i. 

FFTT222( 

FFTT223i 

FFTT224C 

FFTT225; 

FFTT226C 

FFTT227;; 

FFTT228 . 

FFTT229C 

FFTT230C 

FFTT231 : 

FFTT232C 

FFTT233C 

FFTT234C 

FFTT235l 

FFTT236 ; 

FFTT237C 

FFTT238 ; 

FFTT239(; 

FFTT240C 

FFTT241C 

FFTT242C 

FFTT243C 

FFTT244C 

FFTT245C 

FFTT246C 

FFTT247C 

FFTT248C 

FFTT249C 

FFTT250C 

FFTT251C 

FFTT252C 

FFTT2 53f 

FFTT254i; 

FFTT255C 

FFTT256C 

FFTT257: 

FFTT258C 

FFTT259C 

FFTT260C 

FFTT2E1C 

FFTT262C 

FFTT263C 

FFTT264C 

FFTT265 : 

FFTT266C 

FFTT267C 
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DATA( K3) =U1R-U2R 


FFTT263 


DATA( K3+1)=U11-U2I 


FFTT269. 


DATA( K4) =U3R-U4R 


FFTT270> 

520 

DATA( K4+1)=U3I-U4I 


FFTT271 


KMIN = 4=«'( KMIN-J3)+J3 


FFTT272. 


KDIF=KSTEP 


FFTT273. 


IF(KDIF-NP2) 450,530, 530 


FFTT274- 

530 

CONTINUE 


FFTT275i 


M=MMAX-M 


FFTT276( 


IF( ISIGN)540,550,550 


FFTT277 . 

540 

TEMPR=WR 


FFTT278> 


WR=-WI 


FFTT279i 


WI=-TEMPR 


FFTT280. 


GO TO 560 


FFTT281C 

550 

TEMPP=WR 


FFTT282 . 


WR=WI 


FFTT283C 


WI=TEMPR 


FFTT284C 

560 

IF(M-LMAX)565,565,410 


FFTT285:. 

565 

TEMPR=WR 


FFTT286C 


WR = WR*WSTPR-WI=«'WSTPI+WR 


FFTT287C 

570 

WI=WI«WSTPR+TEMPR*WSTPI+WI 


FFTT288 , 


IPAR=3-IPAR 


FFTT289C 


MMAX=MMAX+MMAX 


FFTT290L 


GO TO 360 


FFTT291L 

C 



FFTT292!; 

C 

MAIN LOOP FOR FACTORS NOT EQUAL TO TWO. APPLY THE TWIDDLE 

FACTOR 

FFTT293: 

C 

W = EXP( ISIGN*2*PI*SQRT(-1)*( J2-l)*( J1-J2)/(NP2*IFP1 ) ) , THEN 


FFTT294C 

C 

PERFORM A FOURIER TRANSFORM OF LENGTH IFACT(IF), MAKING USE 

OF 

FFTT295. 

C 

CONJUGATE SYMMETRIES. 


FFTT296 

C 



FFTT297( 

600 

IF (NTWG-NP2 1605,700,700 


FFTT298C 

605 

IFP1=N0N2 


FFTT299 


IF = 1 


FFTT300i 


NPlHF=NPl/2 


FFTT301 . 

610 

IFP2= IFPl/IFACTI IF) 


FFTT302! 


J1RNG=NP2 


FFTT303; 


IF( ICASE-3)612,611,612 


FFTT304. 

611 

J1RNG=(NP2+IFP1 )/2 


FFTT305C 


J2STP = NP2/IFACT( IF) 


FFTT3C6C 


J1RG2=( J2STP+IFP2)/2 


FFTT307. 

612 

J2MIN=1+IFP2 


FFTT303C 


IF( IFP1-NP2) 615,640,640 


FFTT309C 

615 

DO 635 J2=J2MIN, IFPl, IFP2 


FFTTSlOi 


THETA=-TWOPI*FLOAT( J2-1) /FLOAT (NP2) 


FFTT311C 


IF( ISIGN)625,620,620 


FFTT312C 

620 

THETA=-THETA 


FFTT313C 

625 

SINTH=SIN(THETA/2. ) 


FFTT314C 


WSTPR=-2.*SINTH*SINTH 


FFTT315-; 


0 ^^ 
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WSTPI=SIN(THETA) 

WR=WSTPR+1. 

WI=WSTPI 

J1MIN=J2+IFP1 

DO 635 J1 = J1MIN,J1RNG, IFPl 

I lMAX=Jl+IlRNG-2 

DO 630 I1=J1,I1MAX,2 

DO 630 13=11, NT0T,NP2 

J3MAX=I3+IFP2-NP1 

DO 630 J3=I3, J3MAX,NP1 

TEMPR=DATA( J3) 

DATA( J3) =DATA( J3) *WP-OATA( J3+D-WI 
63C DATA( J3+1) =TEMPR*WI+OATA( J3+1)=!'WR 
TEMPR=WR 

WR=WP«WSTPR-WI#WSTP1+WR 
635 WI=TEMPR*WSTPI+WI*WSTPR+WI 
640 THETA=-TWOPI/FLOAT{ I F ACT ( IF) ) 
IF(ISIGN)650,645,645 
645 THETA=-THETA 
650 SlNTH=SIN(THETA/2. ) 

WSTPR=-2 ,«SINTH«SINTH 
WSTPI=SIN(THETA) 

KSTEP=2*N/IFACT( IF) 

KRANG=KSTEP«( IFACT( IF)/2)+l 

00 698 11=1, I1RNG,2 

DO 698 13=11, NT0T,NP2 

DO 690 KMIN=1,KRANG,KSTEP 

J1MAX = I3+J1P.NG-IFP1 

DO 680 J1=I3 , JIMAX, IFPl 

J3MAX=J1+IFP2-NP1 

DO 680 J3=J1, J3MAX,NP1 

J2MAX=J3+IFP1-IFP2 

K= KM I N+( J3- J 1+ I J 1- 1 3 ) / I FACT ( IF ) ) /NP 1 HF 
IF{KMIN-1)655,655,665 
655 SUMR=0. 

SUMI=0. 

DO 660 J2=J3, J2MAX, IFP2 
SUMR=SUMR+DATA(J2) 

660 SUMI=SUMI+DATA(J2+1) 

WORK(K)=SUMR 
WORK! K+1 )=SUMI 
GO TO 680 

665 KC0NJ = K+2’!'(N-KMIN+1 ) 

J2=J2MAX 

SUMR=DATAIJ2) 

SUMI=DATA( J2-H) 

0LDSR=0. 

0LDSI=0. 
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J2=J2-IFP2 



FFTT364 

670 

TEMPR=SUMR 



FFTT365 


TEMPI =SUMI 



FFTT36C 


SUMR = TWOWR=<'SUMR-OLDSR+DATA( J2) 



FFTT367 


SUMI = TWOWR«SUMI-OLDSH-DATA( J2 + 1) 



FFTT366. 


OLDSR=TEMPR 



FFTT369 


OLDSI=TEMPI 



FFTT370;. 


J2=J2-IFP2 



FFTT371; 


IF( J2-J3)675,675,670 



FFTT372 . 

675 

TEMPR=WR#SUMR-0LDSR+0ATAIJ2) 



FFTT373C 


TEMPI=WI*SUMI 



FFTT3741 


WORK( K)=TEMPR-TEMPI 



FFTT375 


W0RK( KCONJ) =TEMPR+TEMPI 



FFTT3 76'.. 


TEMPR=WR*SUMI-OLDSI+OATA( J2+1) 



FFTT377. 


TEMPI=WI*SUMR 



FFT1 3.78. 


W0RK( K+1)=TEMPR+TEMPI 



FFTT379C 


WORK! KCONJ+1 )=TEMPR-TEMPI 



FFTT380 

680 

CONTINUE 



FFTT381i 


IF(KMIN-1)685,685»686 



FFTT382'.. 

685 

WR=WSTPR+1. 



FFTT383. 


WI=WSTPI 



FFTT384C 


GO TO 690 



FFTT385( 

686 

TEMPR=WR 



FFTT386 : 


WR=WR*WSTPR-WI*WSTPI+WR 



FFTT367C 


WI=TEMPR#WSTPI+WI#WSTPR+WI 



FFTT388L 

690 

TWQWR=WR+WR 



FFTT389: 


IF( ICASE-3)692,691»692 



FFTT390C 

691 

IF( IFP1-NP2) 695,692,692 



FFTT391C 

692 

K=1 



FFTT3923 


I2MAX= I3+NP2-NPI 



FFTT393. 


DO 693 12=13, I2MAX,NP1 



FFTT394.. 


DATA { 12) =WORK(K) 



FFTT395C 


DATAI I2+1)=W0RK(K+1) 



FFTT396f 

693 

K=K+2 



FFTT397 


GO TO 698 



FFTT398. 

C 




FFTT399., 

C 

COMPLETE A REAL TRANSFORM IN THE 

1ST DIMENSION, N ODD, 

BY CON- 

FFTT400.. 

C 

JUGATE SYMMETRIES AT EACH STAGE. 



FFTT401' 

C 




FFTT4021 

695 

J3MAX=I3+IFP2-NP1 



FFTT4033 


DO 697 J3=I3, J3MAX,NP1 



■ FFTT404C 


J2MAX=J3+NP2-J2STP 



FFTT405C 


DO 697 J2=J3, J2MAX, J2STP 



FFTT406f. 


J1MAX=J2+J1RG2-IFP2 



FFTT407C 


J1CNJ=J3+J2MAX+J2STP-J2 



FFTT408 ; 


DO 697 J1=J2, JIMAX, IFP2 



FFTT409C 


K=1+J1-I3 



FFTT410C 


DATA( J1)=W0RK(K) 



FFTT411C 
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DATAIJ1+1)=W0RK(K+13 
IF( Jl-J2)697,697,696 

696 DATA( J1CNJ)=W0RK(K) 

DATA( JlCNJ+1 )=-WORK(K+l J 

697 J1CNJ=J1CNJ-IFP2 

698 CONTINUE 
IF=IF+1 
IFP1=IFP2 

IF( IFP1-NP1)700,700,610 
C 

C COMPLETE A REAL TRANSFORM IN THE 1ST DIMENSION^ N EVENt BY CON- 

C JUGATE SYMMETRIES. 

C 

700 GO TO (900,800,900*701) fICASE 

701 NHALF=N 
N=N+N 

THETA=-TWOPI/FLOAT( N) 

IF(ISIGN)703,702,702 

702 THETA=-THETA 

703 SINTH=SIN( THETA/2. ) 

WSTPR=-2.’!'SINTH«SINTH 

WSTPI=SIN(THETA) 

WR=WSTPR + 1. 

WI=WSTPI 

IMIN=3 

JMIN=2*NHALF-1 
GO TO 725 
710 J=JMIN 

DO 720 I=IMIN,NT0T,NP2 
SUMR= (DATA! I )+DATA( J) )/2. 

SUMI =( DATA! I+l )+OATA( J+1) )/2. 

DIFR= (CATA( I )-DATA( J) Ml, 

DIFI=( DATA! I+l )-DATA(J+l) )/2. 

TEMPR=WR*SUMI+WI*DIFR 
TEMPI =WI«SUMI-WR*DIFR 
DATA( I )=SUMR+TEMPR 
DATA( I+l )=DIFI+TEMPI 
DATA( J)=SUMR-TEMPR 
DATA! J+1)=-DIFI+TEMPI 
720 J=J+NP2 

IMIN=IMIN+2 

JMIN=JMIN-2 

TEMPR=WR 

WR=WR«WSTPR-WI*WSTPI+WR 
WI=TEMPR«WSTPI+WI«WSTPR+WI 
725 IF( IMIN-JMIN) 710,730,740 

730 IF(ISIGN)731,740,740 

731 DO 735 I=IMIN,NT0T,NP2 
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735 

740 


745 


750 

755 


760 


765 

770 


775 


730 


C 

C 

c 

c 

800 

805 


810 

820 

830 


840 

850 


OATA( I+l )=-0ATA( I+l) 

NP2=NP2+NP2 

NTOT=NTOT+NTOT 

J=NTOT+l 

IMAX=NTQT/2+l 

IMIN=IMAX-2*NHALF 

I = IMIN 

GO TO 755 

DATAC J)=OATA(I ) 

DATA( J+1)=-DATA(I+1) 

1 = 1 + 2 
J = J-2 

IF( I-IMAX)750,760,760 

DATA( J )=DATA( I M I N ) -OAT A ( I M I N+1 ) 

DATA{ J+1)=0. 

IF(I-J )770,780,780 
DATA( J)=OATA( I ) 

DATA( J + 1J=DATA{ I + l) 

1 = 1-2 
J=J-2 

IF( I-IMIN)775,775,765 

DATA( J ) =DATA { I M I N ) +DAT A ( I MI N+1 ) 

DATA! J+1)=0. 

IMAX= IMIN 
GO TO 745 

DATA( 1)=DATA( 1)+DATA(2) 

0ATA(2)=0. 

GO TO 900 

COMPLETE A REAL TRANSFORM FOP THE 2ND OR 3RD DIMENSION BY 
CONJUGATE SYMMETRIES. 

IF( I1RNG-NP1)805,900,900 

DO 860 I3=1,NT0T,NP2 

I2MAX=I3+NP2-NP1 

DO 860 12=13, I2MAX,NP1 

IMIN=I2+I1RNG 

IMAX=I2+NPl-2 

JMAX=2*I3+NP1-IMIN 

IF( 12-13)820,820,810 

JMAX=JMAX+NP2 

IF( IDIM-2)850,850,830 

J=JMAX+NP0 

DO 840 I=IMIN, IMAX,2 
DATA( I )=DATA{ J) 

DATA( I+1)=-DATA( J+1) 

J = J-2 
J=JMAX 
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DO 860 IMAXtNPO 

OATA( I )=DATA( J) 

DATAi I+li=-DATA{ J+1) 
860 J=J-NPO 


END OF LOOP ON EACH DIMENSION 

900 NP0=NP1 
NP1=NP2 
910 NPREV=N 
920 RETURN 
END 
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